Ignore:
Timestamp:
2020-12-03T10:26:33+01:00 (12 months ago)
Author:
mathiot
Message:

merge tickets_icb_1900 into trunk

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/src/OCE/ICB/icbthm.F90

    r13281 r14030  
    4949      INTEGER, INTENT(in) ::   kt   ! timestep number, just passed to icb_utl_print_berg 
    5050      ! 
    51       INTEGER  ::   ii, ij 
    52       REAL(wp) ::   zM, zT, zW, zL, zSST, zVol, zLn, zWn, zTn, znVol, zIC, zDn 
     51      INTEGER  ::   ii, ij, jk, ikb 
     52      REAL(wp) ::   zM, zT, zW, zL, zSST, zVol, zLn, zWn, zTn, znVol, zIC, zDn, zD, zvb, zub, ztb 
     53      REAL(wp) ::   zMv, zMe, zMb, zmelt, zdvo, zdvob, zdva, zdM, zSs, zdMe, zdMb, zdMv 
    5354      REAL(wp) ::   zSSS, zfzpt 
    54       REAL(wp) ::   zMv, zMe, zMb, zmelt, zdvo, zdva, zdM, zSs, zdMe, zdMb, zdMv 
    5555      REAL(wp) ::   zMnew, zMnew1, zMnew2, zheat_hcflux, zheat_latent, z1_12 
    5656      REAL(wp) ::   zMbits, znMbits, zdMbitsE, zdMbitsM, zLbits, zAbits, zMbb 
    57       REAL(wp) ::   zxi, zyj, zff, z1_rday, z1_e1e2, zdt, z1_dt, z1_dt_e1e2 
     57      REAL(wp) ::   zxi, zyj, zff, z1_rday, z1_e1e2, zdt, z1_dt, z1_dt_e1e2, zdepw 
     58      REAL(wp), DIMENSION(jpk) :: ztoce, zuoce, zvoce, ze3t, zzMv 
    5859      TYPE(iceberg), POINTER ::   this, next 
    5960      TYPE(point)  , POINTER ::   pt 
     
    8586         pt => this%current_point 
    8687         nknberg = this%number(1) 
    87          CALL icb_utl_interp( pt%xi, pt%e1, pt%uo, pt%ui, pt%ua, pt%ssh_x,   & 
    88             &                 pt%yj, pt%e2, pt%vo, pt%vi, pt%va, pt%ssh_y,   & 
    89             &                 pt%sst, pt%cn, pt%hi, zff, pt%sss ) 
     88 
     89         CALL icb_utl_interp( pt%xi, pt%yj,            &   ! position 
     90             &                 pssu=pt%ssu, pua=pt%ua, &   ! oce/atm velocities 
     91             &                 pssv=pt%ssv, pva=pt%va, &   ! oce/atm velocities 
     92             &                 psst=pt%sst, pcn=pt%cn, & 
     93             &                 psss=pt%sss             ) 
     94 
     95         IF ( nn_sample_rate > 0 .AND. MOD(kt-1,nn_sample_rate) == 0 ) THEN 
     96            CALL icb_utl_interp( pt%xi, pt%yj, pe1=pt%e1, pe2=pt%e2,                 & 
     97               &                 pui=pt%ui, pssh_i=pt%ssh_x, & 
     98               &                 pvi=pt%vi, pssh_j=pt%ssh_y, & 
     99               &                 phi=pt%hi,                  & 
     100               &                 plat=pt%lat, plon=pt%lon ) 
     101         END IF 
    90102         ! 
    91103         zSST = pt%sst 
     
    95107         zM   = pt%mass 
    96108         zT   = pt%thickness                               ! total thickness 
    97        ! D   = (rn_rho_bergs/pp_rho_seawater)*zT ! draught (keel depth) 
    98        ! F   = zT - D ! freeboard 
     109         zD   = rho_berg_1_oce * zT                        ! draught (keel depth) 
    99110         zW   = pt%width 
    100111         zL   = pt%length 
     
    108119 
    109120         ! Environment 
    110          zdvo = SQRT( (pt%uvel-pt%uo)**2 + (pt%vvel-pt%vo)**2 ) 
    111          zdva = SQRT( (pt%ua  -pt%uo)**2 + (pt%va  -pt%vo)**2 ) 
    112          zSs  = 1.5_wp * SQRT( zdva ) + 0.1_wp * zdva                ! Sea state      (eqn M.A9) 
    113  
     121         ! default sst, ssu and ssv 
     122         ! ln_M2016: use temp, u and v profile 
     123         IF ( ln_M2016 ) THEN 
     124 
     125            ! load t, u, v and e3 profile at icb position 
     126            CALL icb_utl_interp( pt%xi, pt%yj, ptoce=ztoce, puoce=zuoce, pvoce=zvoce, pe3t=ze3t ) 
     127             
     128            !compute bottom level 
     129            CALL icb_utl_getkb( pt%kb, ze3t, zD ) 
     130 
     131            ikb = MIN(pt%kb,mbkt(ii,ij))                             ! limit pt%kb by mbkt  
     132                                                                     ! => bottom temperature used to fill ztoce(mbkt:jpk) 
     133            ztb = ztoce(ikb)                                         ! basal temperature 
     134            zub = zuoce(ikb) 
     135            zvb = zvoce(ikb) 
     136         ELSE 
     137            ztb = pt%sst 
     138            zub = pt%ssu 
     139            zvb = pt%ssv 
     140         END IF 
     141 
     142         zdvob = SQRT( (pt%uvel-zub)**2 + (pt%vvel-zvb)**2 )        ! relative basal velocity 
     143         zdva  = SQRT( (pt%ua  -pt%ssu)**2 + (pt%va  -pt%ssv)**2 )  ! relative wind 
     144         zSs   = 1.5_wp * SQRT( zdva ) + 0.1_wp * zdva              ! Sea state      (eqn M.A9) 
     145         ! 
    114146         ! Melt rates in m/s (i.e. division by rday) 
    115          zMv = MAX( 7.62d-3*zSST+1.29d-3*(zSST**2)                    , 0._wp ) * z1_rday      ! Buoyant convection at sides (eqn M.A10) 
     147         ! Buoyant convection at sides (eqn M.A10) 
     148         IF ( ln_M2016 ) THEN 
     149            ! averaging along all the iceberg draft 
     150            zzMv(:) = MAX( 7.62d-3*ztoce(:)+1.29d-3*(ztoce(:)**2), 0._wp ) * z1_rday 
     151            CALL icb_utl_zavg(zMv, zzMv, ze3t, zD, ikb ) 
     152         ELSE 
     153            zMv = MAX( 7.62d-3*zSST+1.29d-3*(zSST**2), 0._wp ) * z1_rday 
     154         END IF 
     155         ! 
     156         ! Basal turbulent melting     (eqn M.A7 ) 
    116157         IF ( zSST > zfzpt ) THEN                                                              ! Calculate basal melting only if SST above freezing point   
    117             zMb = MAX( 0.58_wp*(zdvo**0.8_wp)*(zSST+4.0_wp)/(zL**0.2_wp) , 0._wp ) * z1_rday   ! Basal turbulent melting     (eqn M.A7 ) 
     158            zMb = MAX( 0.58_wp*(zdvob**0.8_wp)*(ztb+4.0_wp)/(zL**0.2_wp) , 0._wp ) * z1_rday 
    118159         ELSE 
    119160            zMb = 0._wp                                                                        ! No basal melting if SST below freezing point      
    120161         ENDIF 
    121          zMe = MAX( z1_12*(zSST+2.)*zSs*(1._wp+COS(rpi*(zIC**3)))     , 0._wp ) * z1_rday      ! Wave erosion                (eqn M.A8 ) 
     162         ! 
     163         ! Wave erosion                (eqn M.A8 ) 
     164         zMe = MAX( z1_12*(zSST+2.)*zSs*(1._wp+COS(rpi*(zIC**3)))     , 0._wp ) * z1_rday 
    122165 
    123166         IF( ln_operator_splitting ) THEN      ! Operator split update of volume/mass 
     
    207250 
    208251         ! Rolling 
    209          zDn = ( rn_rho_bergs / pp_rho_seawater ) * zTn       ! draught (keel depth) 
     252         zDn = rho_berg_1_oce * zTn       ! draught (keel depth) 
    210253         IF( zDn > 0._wp .AND. MAX(zWn,zLn) < SQRT( 0.92*(zDn**2) + 58.32*zDn ) ) THEN 
    211254            zT  = zTn 
     
    224267 
    225268!!gm  add a test to avoid over melting ? 
     269!!pm  I agree, over melting could break conservation (more melt than calving) 
    226270 
    227271         IF( zMnew <= 0._wp ) THEN       ! Delete the berg if completely melted 
Note: See TracChangeset for help on using the changeset viewer.