Changeset 13357 for NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbdyn.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/icbdyn.F90 (modified) (5 diffs)

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

@@ -197 +197 @@
       ij  = mj1( ij  )
       !
+      ! assume icb is grounded if tmask(ii,ij,1) or tmask(ii,ij,ikb), depending of the option is not 0
+      !IF ( ln_icb_ground ) THEN
+      !   ! interpol needed data
+      !   CALL icb_utl_interp( pxi, pyj, pe3t=ze3t )   ! 3d velocities
+      !  
+      !   !compute bottom level
+      !   CALL icb_utl_getkb( ikb, ze3t, zD )
+      !
+      !   IF(  tmask(ii,ij,ikb)  /=   0._wp  )   RETURN           ! berg reach a new t-cell, but an ocean one
+      !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
@@ -254 +265 @@
       REAL(wp), PARAMETER ::   pp_Cr0       = 0.06_wp    !
       !
-      INTEGER  ::   itloop
-      REAL(wp) ::   zuo, zui, zua, zuwave, zssh_x, zsst, zcn, zhi
-      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, pyj, pe1=pe1, pe2=pe2,   &   ! scale factor
-         &                 puo=zuo, pui=zui, pua=zua,    &   ! oce/ice/atm velocities
-         &                 pvo=zvo, 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
-
-      IF (lwp) THEN
-         WRITE(numout,*) 'icbdyn ', pxi, pyj 
-         WRITE(numout,*) pe1, zuo, zui, zua, zssh_x
-         WRITE(numout,*) pe2, zvo, zvi, zva, zssh_y
-         WRITE(numout,*) zhi, zff 
-         WRITE(numout,*) ''
-      END IF
-
+      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
@@ -291 +294 @@
       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.
@@ -318 +321 @@
       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 
+         zuo = 0.0 ; zvo = 0.0; zdep = 0.0
+         DO jk = 1, ikb-1
+            zuo = zuo + zuoce(jk) * ze3t(jk)
+            zvo = zvo + zvoce(jk) * ze3t(jk)
+            zdep = zdep + ze3t(jk)
+         END DO
+         zuo = (zuo + zuoce(ikb) * (zD - zdep)) / zD
+         zvo = (zvo + zvoce(ikb) * (zD - zdep)) / zD
+      ELSE
+         zuo = zssu
+         zvo = zssv
+      END IF
+
       zuveln = puvel   ;   zvveln = pvvel ! Copy starting uvel, vvel
       !
@@ -330 +357 @@
          ! 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