Changeset 14360 for NEMO/branches/2021

Timestamp:

2021-01-28T18:15:06+01:00 (3 years ago)

Author:

mathiot

Message:

update speeding ticket code + minor fixes in velocity used in icb_ground for stage 3 and 4

Location:

NEMO/branches/2021/ticket2581_NEMO4.0-HEAD_icb_speeding_ticket

Files:

• cfgs/SHARED/namelist_ref (modified) (1 diff)
• src/OCE/ICB/icbdyn.F90 (modified) (8 diffs)

NEMO/branches/2021/ticket2581_NEMO4.0-HEAD_icb_speeding_ticket/cfgs/SHARED/namelist_ref

@@ -536 +536 @@
    rn_test_box             = 108.0,  116.0, -66.0, -58.0
    ln_use_calving          = .false. ! Use calving data even when nn_test_icebergs > 0
-   rn_speed_limit          = 0.      ! CFL speed limit for a berg
+   rn_speed_limit          = 0.      ! CFL speed limit for a berg (safe value is 0.4, see #2581)
 
    cn_dir      = './'      !  root directory for the calving data location

NEMO/branches/2021/ticket2581_NEMO4.0-HEAD_icb_speeding_ticket/src/OCE/ICB/icbdyn.F90

@@ -14 +14 @@
    USE dom_oce        ! NEMO ocean domain
    USE phycst         ! NEMO physical constants
+   USE in_out_manager                      ! IO parameters
    !
    USE icb_oce        ! define iceberg arrays
@@ -85 +86 @@
          !                                         !**   A1 = A(X1,V1)
          CALL icb_accel( berg , zxi1, ze1, zuvel1, zuvel1, zax1,     &
-            &                   zyj1, ze2, zvvel1, zvvel1, zay1, zdt_2 )
+            &                   zyj1, ze2, zvvel1, zvvel1, zay1, zdt_2, 1./2. )
          !
          zu1 = zuvel1 / ze1                           !**   V1 in d(i,j)/dt
@@ -102 +103 @@
          !                                         !**   A2 = A(X2,V2)
          CALL icb_accel( berg , zxi2, ze1, zuvel2, zuvel1, zax2,    &
-            &                   zyj2, ze2, zvvel2, zvvel1, zay2, zdt_2 )
+            &                   zyj2, ze2, zvvel2, zvvel1, zay2, zdt_2, 1./2. )
          !
          zu2 = zuvel2 / ze1                           !**   V2 in d(i,j)/dt
@@ -113 +114 @@
          zyj3  = zyj1  + zdt_2 * zv2   ;   zvvel3 = zvvel1 + zdt_2 * zay2
          !
-         CALL icb_ground( zxi3, zxi1, zu3,   &
-            &                zyj3, zyj1, zv3, ll_bounced )
+         CALL icb_ground( zxi3, zxi1, zu2,   &
+            &             zyj3, zyj1, zv2, ll_bounced )
 
          !                                         !**   A3 = A(X3,V3)
          CALL icb_accel( berg , zxi3, ze1, zuvel3, zuvel1, zax3,    &
-            &                   zyj3, ze2, zvvel3, zvvel1, zay3, zdt )
+            &                   zyj3, ze2, zvvel3, zvvel1, zay3, zdt, 1. )
          !
          zu3 = zuvel3 / ze1                           !**   V3 in d(i,j)/dt
@@ -129 +130 @@
          zyj4 = zyj1 + zdt * zv3   ;   zvvel4 = zvvel1 + zdt * zay3
 
-         CALL icb_ground( zxi4, zxi1, zu4,   &
-            &             zyj4, zyj1, zv4, ll_bounced )
+         CALL icb_ground( zxi4, zxi1, zu3,   &
+            &             zyj4, zyj1, zv3, ll_bounced )
 
          !                                         !**   A4 = A(X4,V4)
          CALL icb_accel( berg , zxi4, ze1, zuvel4, zuvel1, zax4,    &
-            &                   zyj4, ze2, zvvel4, zvvel1, zay4, zdt )
+            &                   zyj4, ze2, zvvel4, zvvel1, zay4, zdt, 1. )
 
          zu4 = zuvel4 / ze1                           !**   V4 in d(i,j)/dt
@@ -235 +236 @@
 
    SUBROUTINE icb_accel( berg , pxi, pe1, puvel, puvel0, pax,                &
-      &                         pyj, pe2, pvvel, pvvel0, pay, pdt )
+      &                         pyj, pe2, pvvel, pvvel0, pay, pdt, pcfl_scale )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE icb_accel  ***
@@ -244 +245 @@
       !!----------------------------------------------------------------------
       TYPE(iceberg ), POINTER, INTENT(in   ) ::   berg             ! berg
+      REAL(wp)               , INTENT(in   ) ::   pcfl_scale
       REAL(wp)               , INTENT(in   ) ::   pxi   , pyj      ! berg position in (i,j) referential
       REAL(wp)               , INTENT(in   ) ::   puvel , pvvel    ! berg velocity [m/s]
@@ -361 +363 @@
          zspeed = SQRT( zuveln*zuveln + zvveln*zvveln )    ! Speed of berg
          IF( zspeed > 0._wp ) THEN
-            zloc_dx = MIN( pe1, pe2 )                          ! minimum grid spacing
-            zspeed_new = zloc_dx / pdt * rn_speed_limit        ! Speed limit as a factor of dx / dt
+            zloc_dx = MIN( pe1, pe2 )                                ! minimum grid spacing
+            ! cfl scale is function of the RK4 step
+            zspeed_new = zloc_dx / pdt * rn_speed_limit * pcfl_scale ! Speed limit as a factor of dx / dt
             IF( zspeed_new < zspeed ) THEN
-               zuveln = zuveln * ( zspeed_new / zspeed )        ! Scale velocity to reduce speed
-               zvveln = zvveln * ( zspeed_new / zspeed )        ! without changing the direction
+               zuveln = zuveln * ( zspeed_new / zspeed )             ! Scale velocity to reduce speed
+               zvveln = zvveln * ( zspeed_new / zspeed )             ! without changing the direction
+               pax = (zuveln - puvel0)/pdt
+               pay = (zvveln - pvvel0)/pdt
+               WRITE(numout,*) 'speeding ticket (zspeed_new, zspeed): ',zspeed_new, zspeed, pdt, pcfl_scale   
+               CALL FLUSH(numout)
                CALL icb_dia_speed()
             ENDIF