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 14274 – NEMO

Changeset 14274


Ignore:
Timestamp:
2021-01-06T12:23:06+01:00 (3 years ago)
Author:
davestorkey
Message:

UKMO/NEMO_4.0.3_icb_speed_limit2 : science changes.

Location:
NEMO/branches/UKMO/NEMO_4.0.3_icb_speed_limit2/src/OCE/ICB
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/UKMO/NEMO_4.0.3_icb_speed_limit2/src/OCE/ICB/icbdia.F90

    r13587 r14274  
    8383   REAL(wp)                      ::  heat_to_bergs_net, heat_to_ocean_net, melt_net 
    8484   REAL(wp)                      ::  calving_to_bergs_net 
     85   REAL(wp)                      ::  vel_factor_min 
    8586 
    8687   INTEGER                       ::  nbergs_start, nbergs_end, nbergs_calved 
    8788   INTEGER                       ::  nbergs_melted 
    88    INTEGER                       ::  nspeeding_tickets 
     89   INTEGER , DIMENSION(4)        ::  nspeeding_tickets 
    8990   INTEGER , DIMENSION(nclasses) ::  nbergs_calved_by_class 
    9091 
     
    124125      nbergs_calved             = 0 
    125126      nbergs_calved_by_class(:) = 0 
    126       nspeeding_tickets         = 0 
     127      nspeeding_tickets(:)      = 0 
     128      vel_factor_min            = 1._wp 
    127129      stored_heat_end           = 0._wp 
    128130      floating_heat_end         = 0._wp 
     
    157159      IF( lk_mpp ) THEN 
    158160         ALLOCATE( rsumbuf(23) )          ; rsumbuf(:) = 0._wp 
    159          ALLOCATE( nsumbuf(4+nclasses) )  ; nsumbuf(:) = 0 
     161         ALLOCATE( nsumbuf(7+nclasses) )  ; nsumbuf(:) = 0 
    160162         rsumbuf(1) = floating_mass_start 
    161163         rsumbuf(2) = bergs_mass_start 
     
    265267            nsumbuf(2) = nbergs_calved 
    266268            nsumbuf(3) = nbergs_melted 
    267             nsumbuf(4) = nspeeding_tickets 
     269            nsumbuf(4:7) = nspeeding_tickets(:) 
    268270            DO ik = 1, nclasses 
    269                nsumbuf(4+ik) = nbergs_calved_by_class(ik) 
     271               nsumbuf(7+ik) = nbergs_calved_by_class(ik) 
    270272            END DO 
    271             CALL mpp_sum( 'icbdia', nsumbuf(1:nclasses+4), nclasses+4 ) 
     273            CALL mpp_sum( 'icbdia', nsumbuf(1:nclasses+7), nclasses+7 ) 
    272274            ! 
    273275            nbergs_end        = nsumbuf(1) 
    274276            nbergs_calved     = nsumbuf(2) 
    275277            nbergs_melted     = nsumbuf(3) 
    276             nspeeding_tickets = nsumbuf(4) 
     278            nspeeding_tickets(:) = nsumbuf(4:7) 
    277279            DO ik = 1,nclasses 
    278                nbergs_calved_by_class(ik)= nsumbuf(4+ik) 
     280               nbergs_calved_by_class(ik)= nsumbuf(7+ik) 
    279281            END DO 
    280282            ! 
     283            CALL mpp_min( 'icbdia', vel_factor_min, 1 ) 
    281284         ENDIF 
    282285         ! 
     
    329332         IF (nn_verbose_level > 0) THEN 
    330333            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 
     334            WRITE( numicb, '("n speeding tickets by RK4 stage = ",i6,3(",",i6))') (nspeeding_tickets(ik),ik=1,4) 
     335            IF( SUM(nspeeding_tickets) > 0 ) THEN 
     336               WRITE( numicb, '("min velocity reduction factor = ",f12.8)') vel_factor_min 
     337            ENDIF 
    332338         ENDIF 
    333339         ! 
     
    337343         nbergs_calved             = 0 
    338344         nbergs_calved_by_class(:) = 0 
    339          nspeeding_tickets         = 0 
     345         nspeeding_tickets(:)      = 0 
     346         vel_factor_min            = 1._wp 
    340347         stored_heat_start         = stored_heat_end 
    341348         floating_heat_start       = floating_heat_end 
     
    474481 
    475482 
    476    SUBROUTINE icb_dia_speed() 
    477       !!---------------------------------------------------------------------- 
    478       !!---------------------------------------------------------------------- 
    479       ! 
    480       IF( .NOT.ln_bergdia )   RETURN 
    481       nspeeding_tickets = nspeeding_tickets + 1 
     483   SUBROUTINE icb_dia_speed(pvel_factor, pn_stage) 
     484      !!---------------------------------------------------------------------- 
     485      !!---------------------------------------------------------------------- 
     486      REAL(wp), INTENT(in) ::   pvel_factor   ! factor by which velocity reduced 
     487      INTEGER , INTENT(in) ::   pn_stage  ! which stage of the RK4 calculation are we on 
     488      ! 
     489      IF( (.NOT.ln_bergdia) .OR. pn_stage .lt. 1 .OR. pn_stage .gt. 4 )   RETURN 
     490      nspeeding_tickets(pn_stage) = nspeeding_tickets(pn_stage) + 1 
     491      vel_factor_min = MIN(vel_factor_min,pvel_factor)   ! keep track of minimum reduction factor 
    482492      ! 
    483493   END SUBROUTINE icb_dia_speed 
  • NEMO/branches/UKMO/NEMO_4.0.3_icb_speed_limit2/src/OCE/ICB/icbdyn.F90

    r13587 r14274  
    8585         !                                         !**   A1 = A(X1,V1) 
    8686         CALL icb_accel( berg , zxi1, ze1, zuvel1, zuvel1, zax1,     & 
    87             &                   zyj1, ze2, zvvel1, zvvel1, zay1, zdt_2 ) 
     87            &                   zyj1, ze2, zvvel1, zvvel1, zay1, zdt_2, 1 ) 
    8888         ! 
    8989         zu1 = zuvel1 / ze1                           !**   V1 in d(i,j)/dt 
     
    102102         !                                         !**   A2 = A(X2,V2) 
    103103         CALL icb_accel( berg , zxi2, ze1, zuvel2, zuvel1, zax2,    & 
    104             &                   zyj2, ze2, zvvel2, zvvel1, zay2, zdt_2 ) 
     104            &                   zyj2, ze2, zvvel2, zvvel1, zay2, zdt_2, 2 ) 
    105105         ! 
    106106         zu2 = zuvel2 / ze1                           !**   V2 in d(i,j)/dt 
     
    118118         !                                         !**   A3 = A(X3,V3) 
    119119         CALL icb_accel( berg , zxi3, ze1, zuvel3, zuvel1, zax3,    & 
    120             &                   zyj3, ze2, zvvel3, zvvel1, zay3, zdt ) 
     120            &                   zyj3, ze2, zvvel3, zvvel1, zay3, zdt, 3 ) 
    121121         ! 
    122122         zu3 = zuvel3 / ze1                           !**   V3 in d(i,j)/dt 
     
    134134         !                                         !**   A4 = A(X4,V4) 
    135135         CALL icb_accel( berg , zxi4, ze1, zuvel4, zuvel1, zax4,    & 
    136             &                   zyj4, ze2, zvvel4, zvvel1, zay4, zdt ) 
     136            &                   zyj4, ze2, zvvel4, zvvel1, zay4, zdt, 4 ) 
    137137 
    138138         zu4 = zuvel4 / ze1                           !**   V4 in d(i,j)/dt 
     
    235235 
    236236   SUBROUTINE icb_accel( berg , pxi, pe1, puvel, puvel0, pax,                & 
    237       &                         pyj, pe2, pvvel, pvvel0, pay, pdt ) 
     237      &                         pyj, pe2, pvvel, pvvel0, pay, pdt, pstage ) 
    238238      !!---------------------------------------------------------------------- 
    239239      !!                  ***  ROUTINE icb_accel  *** 
     
    250250      REAL(wp)               , INTENT(inout) ::   pax, pay         ! berg acceleration 
    251251      REAL(wp)               , INTENT(in   ) ::   pdt              ! berg time step 
     252      INTEGER                , INTENT(in   ) ::   pstage           ! RK4 stage (for speed limiter) 
    252253      ! 
    253254      REAL(wp), PARAMETER ::   pp_alpha     = 0._wp      ! 
     
    265266      REAL(wp) ::   zampl, zwmod, zCr, zLwavelength, zLcutoff, zLtop 
    266267      REAL(wp) ::   zlambda, zdetA, zA11, zA12, zaxe, zaye, zD_hi 
    267       REAL(wp) ::   zuveln, zvveln, zus, zvs, zspeed, zloc_dx, zspeed_new 
     268      REAL(wp) ::   zuveln, zvveln, zus, zvs, zspeed, zloc_dx, zspeed_new, zvel_factor, zcfl_scale 
    268269      !!---------------------------------------------------------------------- 
    269270 
     
    361362         zspeed = SQRT( zuveln*zuveln + zvveln*zvveln )    ! Speed of berg 
    362363         IF( zspeed > 0._wp ) THEN 
    363             zloc_dx = MIN( pe1, pe2 )                          ! minimum grid spacing 
    364             zspeed_new = zloc_dx / pdt * rn_speed_limit        ! Speed limit as a factor of dx / dt 
     364            zloc_dx = MIN( pe1, pe2 )                                ! minimum grid spacing 
     365            ! cfl scale is function of the RK4 step 
     366            IF( pstage .le. 2 ) THEN 
     367               zcfl_scale=0.5 
     368            ELSE 
     369               zcfl_scale=1.0 
     370            ENDIF 
     371            zspeed_new = zloc_dx / pdt * rn_speed_limit * zcfl_scale ! Speed limit as a factor of dx / dt 
    365372            IF( zspeed_new < zspeed ) THEN 
    366                zuveln = zuveln * ( zspeed_new / zspeed )        ! Scale velocity to reduce speed 
    367                zvveln = zvveln * ( zspeed_new / zspeed )        ! without changing the direction 
    368                CALL icb_dia_speed() 
     373               zvel_factor = zspeed_new / zspeed 
     374               zuveln = zuveln * ( zvel_factor )             ! Scale velocity to reduce speed 
     375               zvveln = zvveln * ( zvel_factor )             ! without changing the direction 
     376               pax = (zuveln - puvel0)/pdt 
     377               pay = (zvveln - pvvel0)/pdt 
     378               WRITE(numicb,*) 'speeding ticket (zspeed_new/zspeed): ',zvel_factor, pdt, zcfl_scale 
     379               CALL FLUSH(numicb) 
     380               CALL icb_dia_speed(zvel_factor,pstage) 
    369381            ENDIF 
    370382         ENDIF 
Note: See TracChangeset for help on using the changeset viewer.