Changeset 14340 for NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/src/OCE/ICB/icbdyn.F90

Timestamp:

2021-01-25T14:59:12+01:00 (3 years ago)

Author:

davestorkey

Message:

Fix iceberg speed limiter. (This is Pierre's version which will probably go back to the NEMO base code at some point).

File:

• NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/src/OCE/ICB/icbdyn.F90 (modified) (8 diffs)

NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/src/OCE/ICB/icbdyn.F90

@@ -85 +85 @@
          !                                         !**   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 )
          !
          zu1 = zuvel1 / ze1                           !**   V1 in d(i,j)/dt
@@ -102 +102 @@
          !                                         !**   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, 2 )
          !
          zu2 = zuvel2 / ze1                           !**   V2 in d(i,j)/dt
@@ -118 +118 @@
          !                                         !**   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, 3 )
          !
          zu3 = zuvel3 / ze1                           !**   V3 in d(i,j)/dt
@@ -134 +134 @@
          !                                         !**   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, 4 )
 
          zu4 = zuvel4 / ze1                           !**   V4 in d(i,j)/dt
@@ -235 +235 @@
 
    SUBROUTINE icb_accel( berg , pxi, pe1, puvel, puvel0, pax,                &
-      &                         pyj, pe2, pvvel, pvvel0, pay, pdt )
+      &                         pyj, pe2, pvvel, pvvel0, pay, pdt, pstage )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE icb_accel  ***
@@ -250 +250 @@
       REAL(wp)               , INTENT(inout) ::   pax, pay         ! berg acceleration
       REAL(wp)               , INTENT(in   ) ::   pdt              ! berg time step
+      INTEGER                , INTENT(in   ) ::   pstage           ! RK4 stage (for speed limiter)
       !
       REAL(wp), PARAMETER ::   pp_alpha     = 0._wp      !
@@ -265 +266 @@
       REAL(wp) ::   zampl, zwmod, zCr, zLwavelength, zLcutoff, zLtop
       REAL(wp) ::   zlambda, zdetA, zA11, zA12, zaxe, zaye, zD_hi
-      REAL(wp) ::   zuveln, zvveln, zus, zvs, zspeed, zloc_dx, zspeed_new
+      REAL(wp) ::   zuveln, zvveln, zus, zvs, zspeed, zloc_dx, zspeed_new, zvel_factor, zcfl_scale
       !!----------------------------------------------------------------------
 
@@ -361 +362 @@
          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
+            IF( pstage .le. 2 ) THEN
+               zcfl_scale=0.5
+            ELSE
+               zcfl_scale=1.0
+            ENDIF
+            zspeed_new = zloc_dx / pdt * rn_speed_limit * zcfl_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
-               CALL icb_dia_speed()
+               zvel_factor = zspeed_new / zspeed
+               zuveln = zuveln * ( zvel_factor )             ! Scale velocity to reduce speed
+               zvveln = zvveln * ( zvel_factor )             ! without changing the direction
+               pax = (zuveln - puvel0)/pdt
+               pay = (zvveln - pvvel0)/pdt
+               WRITE(numicb,*) 'speeding ticket (zspeed_new/zspeed): ',zvel_factor, pdt, zcfl_scale
+               CALL FLUSH(numicb)
+               CALL icb_dia_speed(zvel_factor,pstage)
             ENDIF
          ENDIF