Changeset 13357 for NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbthm.F90

Timestamp:

2020-07-30T13:20:57+02:00 (4 years ago)

Author:

mathiot

Message:

ticket #1900: add changes for Merino 2016 works. Results unchanged if flag ln_M2016 is set to .false. . Remaining step: test ln_M2016 = .true.

File:

• NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbthm.F90 (modified) (5 diffs)

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

@@ -48 +48 @@
       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
-      REAL(wp) ::   zMv, zMe, zMb, zmelt, zdvo, zdva, zdM, zSs, zdMe, zdMb, zdMv
+      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) ::   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, ztmp
       TYPE(iceberg), POINTER ::   this, next
       TYPE(point)  , POINTER ::   pt
@@ -83 +84 @@
          pt => this%current_point
          nknberg = this%number(1)
+
+         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 )
+
          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,                 &
-               &                 puo=pt%uo, pui=pt%ui, pua=pt%ua, pssh_i=pt%ssh_x,   &
-               &                 pvo=pt%vo, pvi=pt%vi, pva=pt%va, pssh_j=pt%ssh_y,   &
-               &                 psst=pt%sst, pcn=pt%cn, phi=pt%hi,                  &
-               &                 plat=pt%lat, plon=pt%lon                            )
-         ELSE
-            CALL icb_utl_interp( pt%xi, pt%yj,           &   ! position
-               &                 puo=pt%uo, pua=pt%ua,   &   ! oce/atm velocities
-               &                 pvo=pt%vo, pva=pt%va,   &   ! oce/atm velocities
-               &                 psst=pt%sst, pcn=pt%cn )    ! sst and ice concentration
-               ! preparation Nacho add t3d and uo, vo, to basal
+               &                 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
          !
@@ -101 +102 @@
          zM   = pt%mass
          zT   = pt%thickness                               ! total thickness
-       ! D   = (rn_rho_bergs/pp_rho_seawater)*zT ! draught (keel depth)
+         zD   = rho_berg_1_oce * zT          ! draught (keel depth)
        ! F   = zT - D ! freeboard
          zW   = pt%width
@@ -114 +115 @@
 
          ! 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))
+            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)
-         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 )
-         zMe = MAX( z1_12*(zSST+2.)*zSs*(1._wp+COS(rpi*(zIC**3)))     , 0._wp ) * z1_rday   ! Wave erosion                (eqn M.A8 )
+         IF ( ln_M2016 ) THEN
+            ! Buoyant convection at sides (eqn M.A10) but averaging along all the iceberg draft
+            ztmp(:) = MAX( 7.62d-3*ztoce(:)+1.29d-3*(ztoce(:)**2), 0._wp ) * z1_rday
+            zMv = 0.0; zdepw = 0.0
+            DO jk = 1,ikb-1
+               zMv = zMv + ze3t(jk)*ztmp(jk)
+               zdepw = zdepw + ze3t(jk)
+            END DO
+            zMv = (zMv + (zD - zdepw)*ztmp(ikb)) / zD
+         ELSE
+            zMv = MAX( 7.62d-3*zSST+1.29d-3*(zSST**2), 0._wp ) * z1_rday
+         END IF
+
+         ! Basal turbulent melting     (eqn M.A7 ) but using the basal temperature
+         zMb = MAX( 0.58_wp*(zdvob**0.8_wp)*(ztb+4.0_wp)/(zL**0.2_wp) , 0._wp ) * z1_rday
+
+         ! 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
@@ -209 +245 @@
 
          ! 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