Ignore:
Timestamp:
2021-02-03T14:12:29+01:00 (6 months ago)
Author:
mathiot
Message:

ticket #2581: implement change made in v4 to trunk version

Location:
NEMO/branches/2021/ticket2581_trunk_icb_speeding_ticket
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2021/ticket2581_trunk_icb_speeding_ticket/cfgs/SHARED/namelist_ref

    r14255 r14382  
    644644   rn_test_box             = 108.0,  116.0, -66.0, -58.0 
    645645   ln_use_calving          = .false. ! Use calving data even when nn_test_icebergs > 0 
    646    rn_speed_limit          = 0.      ! CFL speed limit for a berg 
     646   rn_speed_limit          = 0.      ! CFL speed limit for a berg (safe value is 0.4, see #2581) 
    647647   ! 
    648648   ln_M2016                = .false. ! use Merino et al. (2016) modification (use of 3d ocean data instead of only sea surface data) 
  • NEMO/branches/2021/ticket2581_trunk_icb_speeding_ticket/src/OCE/ICB/icbdia.F90

    r10570 r14382  
    8686   INTEGER                       ::  nbergs_start, nbergs_end, nbergs_calved 
    8787   INTEGER                       ::  nbergs_melted 
    88    INTEGER                       ::  nspeeding_tickets 
     88   INTEGER                       ::  nspeeding_tickets, nspeeding_tickets_all 
    8989   INTEGER , DIMENSION(nclasses) ::  nbergs_calved_by_class 
    9090 
     
    125125      nbergs_calved_by_class(:) = 0 
    126126      nspeeding_tickets         = 0 
     127      nspeeding_tickets_all     = 0 
    127128      stored_heat_end           = 0._wp 
    128129      floating_heat_end         = 0._wp 
     
    271272            CALL mpp_sum( 'icbdia', nsumbuf(1:nclasses+4), nclasses+4 ) 
    272273            ! 
    273             nbergs_end        = nsumbuf(1) 
    274             nbergs_calved     = nsumbuf(2) 
    275             nbergs_melted     = nsumbuf(3) 
    276             nspeeding_tickets = nsumbuf(4) 
     274            nbergs_end            = nsumbuf(1) 
     275            nbergs_calved         = nsumbuf(2) 
     276            nbergs_melted         = nsumbuf(3) 
     277            nspeeding_tickets_all = nsumbuf(4) 
    277278            DO ik = 1,nclasses 
    278279               nbergs_calved_by_class(ik)= nsumbuf(4+ik) 
     
    329330         IF (nn_verbose_level > 0) THEN 
    330331            WRITE( numicb, '("calved by class = ",i6,20(",",i6))') (nbergs_calved_by_class(ik),ik=1,nclasses) 
    331             IF( nspeeding_tickets > 0 )   WRITE( numicb, '("speeding tickets issued = ",i6)') nspeeding_tickets 
     332            IF( nspeeding_tickets_all > 0 ) THEN 
     333                WRITE( numicb, '("speeding tickets issued (this domain)  = ",i6)') nspeeding_tickets 
     334                WRITE( numicb, '("speeding tickets issued (all domains)  = ",i6)') nspeeding_tickets_all 
     335            END IF 
    332336         ENDIF 
    333337         ! 
     
    338342         nbergs_calved_by_class(:) = 0 
    339343         nspeeding_tickets         = 0 
     344         nspeeding_tickets_all     = 0 
    340345         stored_heat_start         = stored_heat_end 
    341346         floating_heat_start       = floating_heat_end 
  • NEMO/branches/2021/ticket2581_trunk_icb_speeding_ticket/src/OCE/ICB/icbdyn.F90

    r14030 r14382  
    8585 
    8686         !                                         !**   A1 = A(X1,V1) 
    87          CALL icb_accel( berg , zxi1, ze1, zuvel1, zuvel1, zax1,     & 
    88             &                   zyj1, ze2, zvvel1, zvvel1, zay1, zdt_2 ) 
     87         CALL icb_accel( kt, berg , zxi1, ze1, zuvel1, zuvel1, zax1,     & 
     88            &                   zyj1, ze2, zvvel1, zvvel1, zay1, zdt_2, 0.5_wp ) 
    8989         ! 
    9090         zu1 = zuvel1 / ze1                           !**   V1 in d(i,j)/dt 
     
    102102 
    103103         !                                         !**   A2 = A(X2,V2) 
    104          CALL icb_accel( berg , zxi2, ze1, zuvel2, zuvel1, zax2,    & 
    105             &                   zyj2, ze2, zvvel2, zvvel1, zay2, zdt_2 ) 
     104         CALL icb_accel( kt, berg , zxi2, ze1, zuvel2, zuvel1, zax2,    & 
     105            &                   zyj2, ze2, zvvel2, zvvel1, zay2, zdt_2, 0.5_wp ) 
    106106         ! 
    107107         zu2 = zuvel2 / ze1                           !**   V2 in d(i,j)/dt 
     
    114114         zyj3  = zyj1  + zdt_2 * zv2   ;   zvvel3 = zvvel1 + zdt_2 * zay2 
    115115         ! 
    116          CALL icb_ground( berg, zxi3, zxi1, zu3,   & 
    117             &                   zyj3, zyj1, zv3, ll_bounced ) 
     116         CALL icb_ground( berg, zxi3, zxi1, zu2,   & 
     117            &                   zyj3, zyj1, zv2, ll_bounced ) 
    118118 
    119119         !                                         !**   A3 = A(X3,V3) 
    120          CALL icb_accel( berg , zxi3, ze1, zuvel3, zuvel1, zax3,    & 
    121             &                   zyj3, ze2, zvvel3, zvvel1, zay3, zdt ) 
     120         CALL icb_accel( kt, berg , zxi3, ze1, zuvel3, zuvel1, zax3,    & 
     121            &                   zyj3, ze2, zvvel3, zvvel1, zay3, zdt, 1._wp ) 
    122122         ! 
    123123         zu3 = zuvel3 / ze1                           !**   V3 in d(i,j)/dt 
     
    130130         zyj4 = zyj1 + zdt * zv3   ;   zvvel4 = zvvel1 + zdt * zay3 
    131131 
    132          CALL icb_ground( berg, zxi4, zxi1, zu4,   & 
    133             &                   zyj4, zyj1, zv4, ll_bounced ) 
     132         CALL icb_ground( berg, zxi4, zxi1, zu3,   & 
     133            &                   zyj4, zyj1, zv3, ll_bounced ) 
    134134 
    135135         !                                         !**   A4 = A(X4,V4) 
    136          CALL icb_accel( berg , zxi4, ze1, zuvel4, zuvel1, zax4,    & 
    137             &                   zyj4, ze2, zvvel4, zvvel1, zay4, zdt ) 
     136         CALL icb_accel( kt, berg , zxi4, ze1, zuvel4, zuvel1, zax4,    & 
     137            &                   zyj4, ze2, zvvel4, zvvel1, zay4, zdt, 1._wp ) 
    138138 
    139139         zu4 = zuvel4 / ze1                           !**   V4 in d(i,j)/dt 
     
    255255 
    256256 
    257    SUBROUTINE icb_accel( berg , pxi, pe1, puvel, puvel0, pax,                & 
    258       &                         pyj, pe2, pvvel, pvvel0, pay, pdt ) 
     257   SUBROUTINE icb_accel( kt, berg , pxi, pe1, puvel, puvel0, pax,                 & 
     258      &                             pyj, pe2, pvvel, pvvel0, pay, pdt, pcfl_scale ) 
    259259      !!---------------------------------------------------------------------- 
    260260      !!                  ***  ROUTINE icb_accel  *** 
     
    265265      !!---------------------------------------------------------------------- 
    266266      TYPE(iceberg ), POINTER, INTENT(in   ) ::   berg             ! berg 
     267      INTEGER                , INTENT(in   ) ::   kt               ! time step 
     268      REAL(wp)               , INTENT(in   ) ::   pcfl_scale 
    267269      REAL(wp)               , INTENT(in   ) ::   pxi   , pyj      ! berg position in (i,j) referential 
    268270      REAL(wp)               , INTENT(in   ) ::   puvel , pvvel    ! berg velocity [m/s] 
     
    404406         zspeed = SQRT( zuveln*zuveln + zvveln*zvveln )    ! Speed of berg 
    405407         IF( zspeed > 0._wp ) THEN 
    406             zloc_dx = MIN( pe1, pe2 )                          ! minimum grid spacing 
    407             zspeed_new = zloc_dx / pdt * rn_speed_limit        ! Speed limit as a factor of dx / dt 
     408            zloc_dx = MIN( pe1, pe2 )                                ! minimum grid spacing 
     409            ! cfl scale is function of the RK4 step 
     410            zspeed_new = zloc_dx / pdt * rn_speed_limit * pcfl_scale ! Speed limit as a factor of dx / dt 
    408411            IF( zspeed_new < zspeed ) THEN 
    409                zuveln = zuveln * ( zspeed_new / zspeed )        ! Scale velocity to reduce speed 
    410                zvveln = zvveln * ( zspeed_new / zspeed )        ! without changing the direction 
     412               zuveln = zuveln * ( zspeed_new / zspeed )             ! Scale velocity to reduce speed 
     413               zvveln = zvveln * ( zspeed_new / zspeed )             ! without changing the direction 
     414               pax = (zuveln - puvel0)/pdt 
     415               pay = (zvveln - pvvel0)/pdt 
     416               ! 
     417               ! print speeding ticket 
     418               IF (nn_verbose_level > 0) THEN 
     419                  WRITE(numicb, 9200) 'icb speeding : ',kt, nknberg, zspeed, & 
     420                       &                pxi, pyj, zuo, zvo, zua, zva, zui, zvi 
     421                  9200 FORMAT(a,i9,i6,f9.2,1x,4(1x,2f9.2)) 
     422               END IF 
     423               ! 
    411424               CALL icb_dia_speed() 
    412425            ENDIF 
  • NEMO/branches/2021/ticket2581_trunk_icb_speeding_ticket/src/OCE/ICB/icbutl.F90

    r14118 r14382  
    320320         ! 
    321321         IF ( ierr > 0 ) THEN 
    322             WRITE(numout,*) 'bottom left corner T point out of bound' 
    323             WRITE(numout,*) pi, kii, mig( 1 ), mig(jpi) 
    324             WRITE(numout,*) pj, kij, mjg( 1 ), mjg(jpj) 
    325             WRITE(numout,*) pmsk 
    326             CALL ctl_stop('STOP','icb_utl_bilin_h: an icebergs coordinates is out of valid range (out of bound error)') 
     322            WRITE(numicb,*) 'bottom left corner T point out of bound' 
     323            WRITE(numicb,*) pi, kii, mig( 1 ), mig(jpi) 
     324            WRITE(numicb,*) pj, kij, mjg( 1 ), mjg(jpj) 
     325            WRITE(numicb,*) pmsk 
     326            CALL FLUSH(numicb) 
     327            CALL ctl_stop('STOP','icb_utl_bilin_e: an icebergs coordinates is out of valid range (out of bound error).'       , & 
     328                 &                                'This can be fixed using rn_speed_limit=0.4 in &namberg.'                   , & 
     329                 &                                'More details in the corresponding iceberg.stat file (nn_verbose_level > 0).' ) 
    327330         END IF 
    328331      END IF 
Note: See TracChangeset for help on using the changeset viewer.