Changeset 14038 for NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/ICB/icbdyn.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/icbdyn.F90 (modified) (14 diffs)

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

@@ -14 +14 @@
    USE dom_oce        ! NEMO ocean domain
    USE phycst         ! NEMO physical constants
+   USE in_out_manager                      ! IO parameters
    !
    USE icb_oce        ! define iceberg arrays
@@ -97 +98 @@
          zyj2 = zyj1 + zdt_2 * zv1          ;   zvvel2 = zvvel1 + zdt_2 * zay1
          !
-         CALL icb_ground( zxi2, zxi1, zu1,   &
-            &             zyj2, zyj1, zv1, ll_bounced )
+         CALL icb_ground( berg, zxi2, zxi1, zu1,   &
+            &                   zyj2, zyj1, zv1, ll_bounced )
 
          !                                         !**   A2 = A(X2,V2)
@@ -113 +114 @@
          zyj3  = zyj1  + zdt_2 * zv2   ;   zvvel3 = zvvel1 + zdt_2 * zay2
          !
-         CALL icb_ground( zxi3, zxi1, zu3,   &
-            &                zyj3, zyj1, zv3, ll_bounced )
+         CALL icb_ground( berg, zxi3, zxi1, zu3,   &
+            &                   zyj3, zyj1, zv3, ll_bounced )
 
          !                                         !**   A3 = A(X3,V3)
@@ -129 +130 @@
          zyj4 = zyj1 + zdt * zv3   ;   zvvel4 = zvvel1 + zdt * zay3
 
-         CALL icb_ground( zxi4, zxi1, zu4,   &
-            &             zyj4, zyj1, zv4, ll_bounced )
+         CALL icb_ground( berg, zxi4, zxi1, zu4,   &
+            &                   zyj4, zyj1, zv4, ll_bounced )
 
          !                                         !**   A4 = A(X4,V4)
@@ -148 +149 @@
          zvvel_n = pt%vvel + zdt_6 * (  zay1 + 2.*(zay2 + zay3) + zay4 )
 
-         CALL icb_ground( zxi_n, zxi1, zuvel_n,   &
-            &             zyj_n, zyj1, zvvel_n, ll_bounced )
+         CALL icb_ground( berg, zxi_n, zxi1, zuvel_n,   &
+            &                   zyj_n, zyj1, zvvel_n, ll_bounced )
 
          pt%uvel = zuvel_n                        !** save in berg structure
@@ -156 +157 @@
          pt%yj   = zyj_n
 
-         ! update actual position
-         pt%lon  = icb_utl_bilin_x(glamt, pt%xi, pt%yj )
-         pt%lat  = icb_utl_bilin(gphit, pt%xi, pt%yj, 'T' )
-
          berg => berg%next                         ! switch to the next berg
          !
@@ -167 +164 @@
 
 
-   SUBROUTINE icb_ground( pi, pi0, pu,   &
-      &                   pj, pj0, pv, ld_bounced )
+   SUBROUTINE icb_ground( berg, pi, pi0, pu,   &
+      &                         pj, pj0, pv, ld_bounced )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE icb_ground  ***
@@ -177 +174 @@
       !!                NB two possibilities available one of which is hard-coded here
       !!----------------------------------------------------------------------
+      TYPE(iceberg ), POINTER, INTENT(in   ) ::   berg             ! berg
+      !
       REAL(wp), INTENT(inout) ::   pi , pj      ! current iceberg position
       REAL(wp), INTENT(in   ) ::   pi0, pj0     ! previous iceberg position
@@ -184 +183 @@
       INTEGER  ::   ii, ii0
       INTEGER  ::   ij, ij0
+      INTEGER  ::   ikb
       INTEGER  ::   ibounce_method
+      !
+      REAL(wp) :: zD 
+      REAL(wp), DIMENSION(jpk) :: ze3t
       !!----------------------------------------------------------------------
       !
@@ -200 +203 @@
       ij  = mj1( ij  )
       !
-      IF(  tmask(ii,ij,1)  /=   0._wp  )   RETURN           ! berg reach a new t-cell, but an ocean one
+      ! assume icb is grounded if tmask(ii,ij,1) or tmask(ii,ij,ikb), depending of the option is not 0
+      IF ( ln_M2016 .AND. ln_icb_grd ) THEN
+         !
+         ! draught (keel depth)
+         zD = rho_berg_1_oce * berg%current_point%thickness
+         !
+         ! interpol needed data
+         CALL icb_utl_interp( pi, pj, pe3t=ze3t )
+         ! 
+         !compute bottom level
+         CALL icb_utl_getkb( ikb, ze3t, zD )
+         !
+         ! berg reach a new t-cell, but an ocean one
+         ! .AND. needed in case berg hit an isf (tmask(ii,ij,1) == 0 and tmask(ii,ij,ikb) /= 0)
+         IF(  tmask(ii,ij,ikb) /= 0._wp .AND. tmask(ii,ij,1) /= 0._wp ) RETURN
+         !
+      ELSE
+         IF(  tmask(ii,ij,1)  /=   0._wp  )   RETURN           ! berg reach a new t-cell, but an ocean one
+      END IF
       !
       ! From here, berg have reach land: treat grounding/bouncing
@@ -257 +278 @@
       REAL(wp), PARAMETER ::   pp_Cr0       = 0.06_wp    !
       !
-      INTEGER  ::   itloop
-      REAL(wp) ::   zuo, zui, zua, zuwave, zssh_x, zsst, zcn, zhi, zsss
-      REAL(wp) ::   zvo, zvi, zva, zvwave, zssh_y
+      INTEGER  ::   itloop, ikb, jk
+      REAL(wp) ::   zuo, zssu, zui, zua, zuwave, zssh_x, zcn, zhi
+      REAL(wp) ::   zvo, zssv, zvi, zva, zvwave, zssh_y
       REAL(wp) ::   zff, zT, zD, zW, zL, zM, zF
       REAL(wp) ::   zdrag_ocn, zdrag_atm, zdrag_ice, zwave_rad
-      REAL(wp) ::   z_ocn, z_atm, z_ice
+      REAL(wp) ::   z_ocn, z_atm, z_ice, zdep
       REAL(wp) ::   zampl, zwmod, zCr, zLwavelength, zLcutoff, zLtop
       REAL(wp) ::   zlambda, zdetA, zA11, zA12, zaxe, zaye, zD_hi
       REAL(wp) ::   zuveln, zvveln, zus, zvs, zspeed, zloc_dx, zspeed_new
+      REAL(wp), DIMENSION(jpk) :: zuoce, zvoce, ze3t, zdepw
       !!----------------------------------------------------------------------
 
       ! Interpolate gridded fields to berg
       nknberg = berg%number(1)
-      CALL icb_utl_interp( pxi, pe1, zuo, zui, zua, zssh_x,                     &
-         &                 pyj, pe2, zvo, zvi, zva, zssh_y, zsst, zcn, zhi, zff, zsss )
+      CALL icb_utl_interp( pxi, pyj, pe1=pe1, pe2=pe2,     &   ! scale factor
+         &                 pssu=zssu, pui=zui, pua=zua,    &   ! oce/ice/atm velocities
+         &                 pssv=zssv, pvi=zvi, pva=zva,    &   ! oce/ice/atm velocities
+         &                 pssh_i=zssh_x, pssh_j=zssh_y,   &   ! ssh gradient
+         &                 phi=zhi, pff=zff)                   ! ice thickness and coriolis
 
       zM = berg%current_point%mass
       zT = berg%current_point%thickness               ! total thickness
-      zD = ( rn_rho_bergs / pp_rho_seawater ) * zT    ! draught (keel depth)
+      zD = rho_berg_1_oce * zT                        ! draught (keel depth)
       zF = zT - zD                                    ! freeboard
       zW = berg%current_point%width
@@ -282 +307 @@
       zhi   = MIN( zhi   , zD    )
       zD_hi = MAX( 0._wp, zD-zhi )
-
-      ! Wave radiation
-      zuwave = zua - zuo   ;   zvwave = zva - zvo     ! Use wind speed rel. to ocean for wave model
+ 
+     ! Wave radiation
+      zuwave = zua - zssu   ;   zvwave = zva - zssv   ! Use wind speed rel. to ocean for wave model
       zwmod  = zuwave*zuwave + zvwave*zvwave          ! The wave amplitude and length depend on the  current;
       !                                               ! wind speed relative to the ocean. Actually wmod is wmod**2 here.
@@ -309 +334 @@
       IF( abs(zui) + abs(zvi) == 0._wp )   z_ice = 0._wp
 
+      ! lateral velocities
+      ! default ssu and ssv
+      ! ln_M2016: mean velocity along the profile
+      IF ( ln_M2016 ) THEN
+         ! interpol needed data
+         CALL icb_utl_interp( pxi, pyj, puoce=zuoce, pvoce=zvoce, pe3t=ze3t )   ! 3d velocities
+        
+         !compute bottom level
+         CALL icb_utl_getkb( ikb, ze3t, zD )
+         
+         ! compute mean velocity 
+         CALL icb_utl_zavg(zuo, zuoce, ze3t, zD, ikb)
+         CALL icb_utl_zavg(zvo, zvoce, ze3t, zD, ikb)
+      ELSE
+         zuo = zssu
+         zvo = zssv
+      END IF
+
       zuveln = puvel   ;   zvveln = pvvel ! Copy starting uvel, vvel
       !
@@ -321 +364 @@
          ! Explicit accelerations
          !zaxe= zff*pvvel -grav*zssh_x +zwave_rad*zuwave &
-         !    -zdrag_ocn*(puvel-zuo) -zdrag_atm*(puvel-zua) -zdrag_ice*(puvel-zui)
+         !    -zdrag_ocn*(puvel-zssu) -zdrag_atm*(puvel-zua) -zdrag_ice*(puvel-zui)
          !zaye=-zff*puvel -grav*zssh_y +zwave_rad*zvwave &
-         !    -zdrag_ocn*(pvvel-zvo) -zdrag_atm*(pvvel-zva) -zdrag_ice*(pvvel-zvi)
+         !    -zdrag_ocn*(pvvel-zssv) -zdrag_atm*(pvvel-zva) -zdrag_ice*(pvvel-zvi)
          zaxe = -grav * zssh_x + zwave_rad * zuwave
          zaye = -grav * zssh_y + zwave_rad * zvwave