Changeset 14038 for NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/ICB/icbthm.F90

Timestamp:

2020-12-03T12:25:19+01:00 (3 years ago)

Author:

laurent

Message:

Catch up with trunk at rev r14037

File:

• NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/ICB/icbthm.F90 (modified) (6 diffs)

NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/ICB/icbthm.F90

@@ -49 +49 @@
       INTEGER, INTENT(in) ::   kt   ! timestep number, just passed to icb_utl_print_berg
       !
-      INTEGER  ::   ii, ij
-      REAL(wp) ::   zM, zT, zW, zL, zSST, zVol, zLn, zWn, zTn, znVol, zIC, zDn
+      INTEGER  ::   ii, ij, jk, ikb
+      REAL(wp) ::   zM, zT, zW, zL, zSST, zVol, zLn, zWn, zTn, znVol, zIC, zDn, zD, zvb, zub, ztb
+      REAL(wp) ::   zMv, zMe, zMb, zmelt, zdvo, zdvob, zdva, zdM, zSs, zdMe, zdMb, zdMv
       REAL(wp) ::   zSSS, zfzpt
-      REAL(wp) ::   zMv, zMe, zMb, zmelt, zdvo, zdva, zdM, zSs, zdMe, zdMb, zdMv
       REAL(wp) ::   zMnew, zMnew1, zMnew2, zheat_hcflux, zheat_latent, z1_12
       REAL(wp) ::   zMbits, znMbits, zdMbitsE, zdMbitsM, zLbits, zAbits, zMbb
-      REAL(wp) ::   zxi, zyj, zff, z1_rday, z1_e1e2, zdt, z1_dt, z1_dt_e1e2
+      REAL(wp) ::   zxi, zyj, zff, z1_rday, z1_e1e2, zdt, z1_dt, z1_dt_e1e2, zdepw
+      REAL(wp), DIMENSION(jpk) :: ztoce, zuoce, zvoce, ze3t, zzMv
       TYPE(iceberg), POINTER ::   this, next
       TYPE(point)  , POINTER ::   pt
@@ -85 +86 @@
          pt => this%current_point
          nknberg = this%number(1)
-         CALL icb_utl_interp( pt%xi, pt%e1, pt%uo, pt%ui, pt%ua, pt%ssh_x,   &
-            &                 pt%yj, pt%e2, pt%vo, pt%vi, pt%va, pt%ssh_y,   &
-            &                 pt%sst, pt%cn, pt%hi, zff, pt%sss )
+
+         CALL icb_utl_interp( pt%xi, pt%yj,            &   ! position
+             &                 pssu=pt%ssu, pua=pt%ua, &   ! oce/atm velocities
+             &                 pssv=pt%ssv, pva=pt%va, &   ! oce/atm velocities
+             &                 psst=pt%sst, pcn=pt%cn, &
+             &                 psss=pt%sss             )
+
+         IF ( nn_sample_rate > 0 .AND. MOD(kt-1,nn_sample_rate) == 0 ) THEN
+            CALL icb_utl_interp( pt%xi, pt%yj, pe1=pt%e1, pe2=pt%e2,                 &
+               &                 pui=pt%ui, pssh_i=pt%ssh_x, &
+               &                 pvi=pt%vi, pssh_j=pt%ssh_y, &
+               &                 phi=pt%hi,                  &
+               &                 plat=pt%lat, plon=pt%lon )
+         END IF
          !
          zSST = pt%sst
@@ -95 +107 @@
          zM   = pt%mass
          zT   = pt%thickness                               ! total thickness
-       ! D   = (rn_rho_bergs/pp_rho_seawater)*zT ! draught (keel depth)
-       ! F   = zT - D ! freeboard
+         zD   = rho_berg_1_oce * zT                        ! draught (keel depth)
          zW   = pt%width
          zL   = pt%length
@@ -108 +119 @@
 
          ! Environment
-         zdvo = SQRT( (pt%uvel-pt%uo)**2 + (pt%vvel-pt%vo)**2 )
-         zdva = SQRT( (pt%ua  -pt%uo)**2 + (pt%va  -pt%vo)**2 )
-         zSs  = 1.5_wp * SQRT( zdva ) + 0.1_wp * zdva                ! Sea state      (eqn M.A9)
-
+         ! default sst, ssu and ssv
+         ! ln_M2016: use temp, u and v profile
+         IF ( ln_M2016 ) THEN
+
+            ! load t, u, v and e3 profile at icb position
+            CALL icb_utl_interp( pt%xi, pt%yj, ptoce=ztoce, puoce=zuoce, pvoce=zvoce, pe3t=ze3t )
+            
+            !compute bottom level
+            CALL icb_utl_getkb( pt%kb, ze3t, zD )
+
+            ikb = MIN(pt%kb,mbkt(ii,ij))                             ! limit pt%kb by mbkt 
+                                                                     ! => bottom temperature used to fill ztoce(mbkt:jpk)
+            ztb = ztoce(ikb)                                         ! basal temperature
+            zub = zuoce(ikb)
+            zvb = zvoce(ikb)
+         ELSE
+            ztb = pt%sst
+            zub = pt%ssu
+            zvb = pt%ssv
+         END IF
+
+         zdvob = SQRT( (pt%uvel-zub)**2 + (pt%vvel-zvb)**2 )        ! relative basal velocity
+         zdva  = SQRT( (pt%ua  -pt%ssu)**2 + (pt%va  -pt%ssv)**2 )  ! relative wind
+         zSs   = 1.5_wp * SQRT( zdva ) + 0.1_wp * zdva              ! Sea state      (eqn M.A9)
+         !
          ! Melt rates in m/s (i.e. division by rday)
-         zMv = MAX( 7.62d-3*zSST+1.29d-3*(zSST**2)                    , 0._wp ) * z1_rday      ! Buoyant convection at sides (eqn M.A10)
+         ! Buoyant convection at sides (eqn M.A10)
+         IF ( ln_M2016 ) THEN
+            ! averaging along all the iceberg draft
+            zzMv(:) = MAX( 7.62d-3*ztoce(:)+1.29d-3*(ztoce(:)**2), 0._wp ) * z1_rday
+            CALL icb_utl_zavg(zMv, zzMv, ze3t, zD, ikb )
+         ELSE
+            zMv = MAX( 7.62d-3*zSST+1.29d-3*(zSST**2), 0._wp ) * z1_rday
+         END IF
+         !
+         ! Basal turbulent melting     (eqn M.A7 )
          IF ( zSST > zfzpt ) THEN                                                              ! Calculate basal melting only if SST above freezing point  
-            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 )
+            zMb = MAX( 0.58_wp*(zdvob**0.8_wp)*(ztb+4.0_wp)/(zL**0.2_wp) , 0._wp ) * z1_rday
          ELSE
             zMb = 0._wp                                                                        ! No basal melting if SST below freezing point     
          ENDIF
-         zMe = MAX( z1_12*(zSST+2.)*zSs*(1._wp+COS(rpi*(zIC**3)))     , 0._wp ) * z1_rday      ! Wave erosion                (eqn M.A8 )
+         !
+         ! Wave erosion                (eqn M.A8 )
+         zMe = MAX( z1_12*(zSST+2.)*zSs*(1._wp+COS(rpi*(zIC**3)))     , 0._wp ) * z1_rday
 
          IF( ln_operator_splitting ) THEN      ! Operator split update of volume/mass
@@ -207 +250 @@
 
          ! Rolling
-         zDn = ( rn_rho_bergs / pp_rho_seawater ) * zTn       ! draught (keel depth)
+         zDn = rho_berg_1_oce * zTn       ! draught (keel depth)
          IF( zDn > 0._wp .AND. MAX(zWn,zLn) < SQRT( 0.92*(zDn**2) + 58.32*zDn ) ) THEN
             zT  = zTn
@@ -224 +267 @@
 
 !!gm  add a test to avoid over melting ?
+!!pm  I agree, over melting could break conservation (more melt than calving)
 
          IF( zMnew <= 0._wp ) THEN       ! Delete the berg if completely melted