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 13359 for NEMO – NEMO

Changeset 13359 for NEMO


Ignore:
Timestamp:
2020-07-30T15:40:36+02:00 (4 years ago)
Author:
mathiot
Message:

ticket #1900: add subroutine icb_utl_zavg, fix issue in icb_utl_getkb, add subroutine test_icb_utl_getkb to test awkward cases (will be removed after the review), simplify icbdyn.F90 icbthm.F90 according to the new routine created.

Location:
NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbdyn.F90

    r13357 r13359  
    332332          
    333333         ! compute mean velocity  
    334          zuo = 0.0 ; zvo = 0.0; zdep = 0.0 
    335          DO jk = 1, ikb-1 
    336             zuo = zuo + zuoce(jk) * ze3t(jk) 
    337             zvo = zvo + zvoce(jk) * ze3t(jk) 
    338             zdep = zdep + ze3t(jk) 
    339          END DO 
    340          zuo = (zuo + zuoce(ikb) * (zD - zdep)) / zD 
    341          zvo = (zvo + zvoce(ikb) * (zD - zdep)) / zD 
     334         CALL icb_utl_zavg(zuo, zuoce, ze3t, zD, ikb) 
     335         CALL icb_utl_zavg(zvo, zvoce, ze3t, zD, ikb) 
    342336      ELSE 
    343337         zuo = zssu 

  • NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbstp.F90

    r13357 r13359  
    7171      ! 
    7272      nktberg = kt 
     73      ! 
     74      !CALL test_icb_utl_getkb 
     75      !CALL ctl_stop('end test icb') 
    7376      ! 
    7477      IF( nn_test_icebergs < 0 .OR. ln_use_calving ) THEN !* read calving data 

  • NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbthm.F90

    r13357 r13359  
    5454      REAL(wp) ::   zMbits, znMbits, zdMbitsE, zdMbitsM, zLbits, zAbits, zMbb 
    5555      REAL(wp) ::   zxi, zyj, zff, z1_rday, z1_e1e2, zdt, z1_dt, z1_dt_e1e2, zdepw 
    56       REAL(wp), DIMENSION(jpk) :: ztoce, zuoce, zvoce, ze3t, ztmp 
     56      REAL(wp), DIMENSION(jpk) :: ztoce, zuoce, zvoce, ze3t, zzMv 
    5757      TYPE(iceberg), POINTER ::   this, next 
    5858      TYPE(point)  , POINTER ::   pt 
    … …  
    142142         IF ( ln_M2016 ) THEN 
    143143            ! Buoyant convection at sides (eqn M.A10) but averaging along all the iceberg draft 
    144             ztmp(:) = MAX( 7.62d-3*ztoce(:)+1.29d-3*(ztoce(:)**2), 0._wp ) * z1_rday 
    145             zMv = 0.0; zdepw = 0.0 
    146             DO jk = 1,ikb-1 
    147                zMv = zMv + ze3t(jk)*ztmp(jk) 
    148                zdepw = zdepw + ze3t(jk) 
    149             END DO 
    150             zMv = (zMv + (zD - zdepw)*ztmp(ikb)) / zD 
     144            zzMv(:) = MAX( 7.62d-3*ztoce(:)+1.29d-3*(ztoce(:)**2), 0._wp ) * z1_rday 
     145            CALL icb_utl_zavg(zMv, zzMv, ze3t, zD, ikb ) 
    151146         ELSE 
    152147            zMv = MAX( 7.62d-3*zSST+1.29d-3*(zSST**2), 0._wp ) * z1_rday 
    … …  
    262257 
    263258!!gm  add a test to avoid over melting ? 
     259!!pm  I agree, over melting could break conservation (more melt than calving) 
    264260 
    265261         IF( zMnew <= 0._wp ) THEN       ! Delete the berg if completely melted 

  • NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbutl.F90

    r13357 r13359  
    4141   PUBLIC   icb_utl_copy          ! routine called in icbstp module 
    4242   PUBLIC   icb_utl_getkb         ! routine called in icbdyn and icbthm modules 
     43   PUBLIC   test_icb_utl_getkb    ! routine called in icbdyn and icbthm modules 
     44   PUBLIC   icb_utl_zavg          ! routine called in icbdyn and icbthm modules 
    4345   PUBLIC   icb_utl_interp        ! routine called in icbdyn, icbthm modules 
    4446   PUBLIC   icb_utl_bilin_h       ! routine called in icbdyn module 
    … …  
    482484      !! 
    483485      !!---------------------------------------------------------------------- 
    484       INTEGER, INTENT(out)               :: kb 
     486      INTEGER,                INTENT(out):: kb 
    485487      REAL(wp), DIMENSION(:), INTENT(in) :: pe3 
    486488      REAL(wp),               INTENT(in) :: pD 
    … …  
    488490      INTEGER  :: jk 
    489491      REAL(wp) :: zdepw 
    490       !! 
    491       zdepw = 0.0 
    492       kb = 1 
     492      !!---------------------------------------------------------------------- 
     493      !! 
     494      zdepw = pe3(1) ; kb = 2 
    493495      DO WHILE ( zdepw <  pD) 
    494496         zdepw = zdepw + pe3(kb) 
    495497         kb = kb + 1 
    496498      END DO 
    497       kb = kb - 1 
     499      kb = MIN(kb - 1,jpk) 
     500   END SUBROUTINE 
     501 
     502   SUBROUTINE icb_utl_zavg(pzavg, pdat, pe3, pD, kb ) 
     503      !!---------------------------------------------------------------------- 
     504      !!                ***  ROUTINE icb_utl_getkb         *** 
     505      !! 
     506      !! ** Purpose :   compute the vertical average of ocean properties affected by icb 
     507      !! 
     508      !!---------------------------------------------------------------------- 
     509      INTEGER,                INTENT(in ) :: kb        ! deepest level affected by icb 
     510      REAL(wp), DIMENSION(:), INTENT(in ) :: pe3, pdat ! vertical profile 
     511      REAL(wp),               INTENT(in ) :: pD        ! draft 
     512      REAL(wp),               INTENT(out) :: pzavg     ! z average 
     513      !!---------------------------------------------------------------------- 
     514      INTEGER  :: jk 
     515      REAL(wp) :: zdep 
     516      !!---------------------------------------------------------------------- 
     517      pzavg = 0.0 ; zdep = 0.0 
     518      DO jk = 1,kb-1 
     519         pzavg = pzavg + pe3(jk)*pdat(jk) 
     520         zdep  = zdep  + pe3(jk) 
     521      END DO 
     522      ! if kb is limited by mbkt  => bottom value is used between bathy and icb tail 
     523      ! if kb not limited by mbkt => ocean value over mask is used (ie 0.0 for u, v) 
     524      pzavg = ( pzavg + (pD - zdep)*pdat(kb)) / pD 
    498525   END SUBROUTINE 
    499526 
    … …  
    855882   END FUNCTION icb_utl_heat 
    856883 
     884   SUBROUTINE test_icb_utl_getkb 
     885      INTEGER :: ikb 
     886      REAL(wp) :: zD, zout 
     887      REAL(wp), DIMENSION(jpk) :: ze3, zin 
     888      WRITE(numout,*) 'Test icb_utl_getkb : ' 
     889      zD = 0.0 ; ze3= 20.0 
     890      WRITE(numout,*) 'INPUT : zD = ',zD,' ze3 = ',ze3(1) 
     891      CALL icb_utl_getkb(ikb, ze3, zD) 
     892      WRITE(numout,*) 'OUTPUT : kb = ',ikb 
     893 
     894      zD = 8000000.0 ; ze3= 20.0 
     895      WRITE(numout,*) 'INPUT : zD = ',zD,' ze3 = ',ze3(1) 
     896      CALL icb_utl_getkb(ikb, ze3, zD) 
     897      WRITE(numout,*) 'OUTPUT : kb = ',ikb 
     898 
     899      zD = 80.0 ; ze3= 20.0 
     900      WRITE(numout,*) 'INPUT : zD = ',zD,' ze3 = ',ze3(1) 
     901      CALL icb_utl_getkb(ikb, ze3, zD) 
     902      WRITE(numout,*) 'OUTPUT : kb = ',ikb 
     903 
     904      zD = 85.0 ; ze3= 20.0 
     905      WRITE(numout,*) 'INPUT : zD = ',zD,' ze3 = ',ze3(1) 
     906      CALL icb_utl_getkb(ikb, ze3, zD) 
     907      WRITE(numout,*) 'OUTPUT : kb = ',ikb 
     908 
     909      zD = 75.0 ; ze3= 20.0 
     910      WRITE(numout,*) 'INPUT : zD = ',zD,' ze3 = ',ze3(1) 
     911      CALL icb_utl_getkb(ikb, ze3, zD) 
     912      WRITE(numout,*) 'OUTPUT : kb = ',ikb 
     913 
     914      WRITE(numout,*) '==================================' 
     915      WRITE(numout,*) 'Test icb_utl_zavg' 
     916      zD = 0.0 ; ze3= 20.0 ; zin=1.0 
     917      CALL icb_utl_getkb(ikb, ze3, zD) 
     918      CALL icb_utl_zavg(zout, zin, ze3, zD, ikb) 
     919      WRITE(numout,*) 'INPUT  : zD = ',zD,' ze3 = ',ze3(1),' zin = ', zin, ' ikb = ',ikb 
     920      WRITE(numout,*) 'OUTPUT : zout = ',zout 
     921 
     922      zD = 50.0 ; ze3= 20.0 ; zin=1.0; zin(3:jpk) = 0.0 
     923      CALL icb_utl_getkb(ikb, ze3, zD) 
     924      CALL icb_utl_zavg(zout, zin, ze3, zD, ikb) 
     925      WRITE(numout,*) 'INPUT  : zD = ',zD,' ze3 = ',ze3(1),' zin = ', zin, ' ikb = ',ikb 
     926      WRITE(numout,*) 'OUTPUT : zout = ',zout 
     927      CALL FLUSH(numout) 
     928 
     929      zD = 80.0 ; ze3= 20.0 ; zin=1.0; zin(3:jpk) = 0.0 
     930      CALL icb_utl_getkb(ikb, ze3, zD) 
     931      CALL icb_utl_zavg(zout, zin, ze3, zD, ikb) 
     932      WRITE(numout,*) 'INPUT  : zD = ',zD,' ze3 = ',ze3(1),' zin = ', zin, ' ikb = ',ikb 
     933      WRITE(numout,*) 'OUTPUT : zout = ',zout 
     934 
     935      zD = 80 ; ze3= 20.0 ; zin=1.0 ; zin(3:jpk) = 0.0 
     936      CALL icb_utl_getkb(ikb, ze3, zD) 
     937      ikb = 2 
     938      CALL icb_utl_zavg(zout, zin, ze3, zD, ikb) 
     939      WRITE(numout,*) 'INPUT  : zD = ',zD,' ze3 = ',ze3(1),' zin = ', zin, ' ikb = ',ikb 
     940      WRITE(numout,*) 'OUTPUT : zout = ',zout 
     941 
     942      CALL FLUSH(numout) 
     943 
     944   END SUBROUTINE test_icb_utl_getkb 
     945 
    857946   !!====================================================================== 
    858947END MODULE icbutl 
Note: See TracChangeset for help on using the changeset viewer.