Changeset 11897

Timestamp:

2019-11-13T13:27:19+01:00 (5 years ago)

Author:

cguiavarch

Message:

Change to add basal melt only if the SST is above the freezing point

Location:

NEMO/branches/UKMO/NEMO_4.0.1_ICB_melting_temperature/src/OCE/ICB

Files:

• icb_oce.F90 (modified) (3 diffs)
• icbdyn.F90 (modified) (2 diffs)
• icbini.F90 (modified) (1 diff)
• icbthm.F90 (modified) (4 diffs)
• icbutl.F90 (modified) (5 diffs)

NEMO/branches/UKMO/NEMO_4.0.1_ICB_melting_temperature/src/OCE/ICB/icb_oce.F90

@@ -57 +57 @@
    TYPE, PUBLIC ::   point              !: properties of an individual iceberg (position, mass, size, etc...)
       INTEGER  ::   year
-      REAL(wp) ::   xi , yj                                              ! iceberg coordinates in the (i,j) referential (global)
-      REAL(wp) ::   e1 , e2                                              ! horizontal scale factors at the iceberg position
-      REAL(wp) ::   lon, lat, day                                        ! geographic position
-      REAL(wp) ::   mass, thickness, width, length, uvel, vvel           ! iceberg physical properties
-      REAL(wp) ::   uo, vo, ui, vi, ua, va, ssh_x, ssh_y, sst, cn, hi    ! properties of iceberg environment 
+      REAL(wp) ::   xi , yj                                                   ! iceberg coordinates in the (i,j) referential (global)
+      REAL(wp) ::   e1 , e2                                                   ! horizontal scale factors at the iceberg position
+      REAL(wp) ::   lon, lat, day                                             ! geographic position
+      REAL(wp) ::   mass, thickness, width, length, uvel, vvel                ! iceberg physical properties
+      REAL(wp) ::   uo, vo, ui, vi, ua, va, ssh_x, ssh_y, sst, cn, hi, sss    ! properties of iceberg environment 
       REAL(wp) ::   mass_of_bits, heat_density
    END TYPE point
@@ -86 +86 @@
    ! particularly for MPP when iceberg can lie inside T grid but outside U, V, or f grid
    REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   uo_e, vo_e
-   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ff_e, tt_e, fr_e
+   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ff_e, tt_e, fr_e, ss_e
    REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ua_e, va_e
    REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ssh_e
@@ -179 +179 @@
          &      ff_e(0:jpi+1,0:jpj+1) , fr_e(0:jpi+1,0:jpj+1)  ,   &
          &      tt_e(0:jpi+1,0:jpj+1) , ssh_e(0:jpi+1,0:jpj+1) ,   &
+         &      ss_e(0:jpi+1,0:jpj+1) ,                            & 
          &      first_width(nclasses) , first_length(nclasses) ,   &
          &      src_calving (jpi,jpj) ,                            &

NEMO/branches/UKMO/NEMO_4.0.1_ICB_melting_temperature/src/OCE/ICB/icbdyn.F90

@@ -258 +258 @@
       !
       INTEGER  ::   itloop
-      REAL(wp) ::   zuo, zui, zua, zuwave, zssh_x, zsst, zcn, zhi
+      REAL(wp) ::   zuo, zui, zua, zuwave, zssh_x, zsst, zcn, zhi, zsss
       REAL(wp) ::   zvo, zvi, zva, zvwave, zssh_y
       REAL(wp) ::   zff, zT, zD, zW, zL, zM, zF
@@ -271 +271 @@
       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 )
+         &                 pyj, pe2, zvo, zvi, zva, zssh_y, zsst, zcn, zhi, zff, zsss )
 
       zM = berg%current_point%mass

NEMO/branches/UKMO/NEMO_4.0.1_ICB_melting_temperature/src/OCE/ICB/icbini.F90

@@ -79 +79 @@
       ua_e(:,:) = 0._wp   ;   va_e(:,:) = 0._wp   ;
       ff_e(:,:) = 0._wp   ;   tt_e(:,:) = 0._wp   ;
-      fr_e(:,:) = 0._wp   ;
+      fr_e(:,:) = 0._wp   ;   ss_e(:,:) = 0._wp   ;
 #if defined key_si3
       hi_e(:,:) = 0._wp   ;

NEMO/branches/UKMO/NEMO_4.0.1_ICB_melting_temperature/src/OCE/ICB/icbthm.F90

@@ -20 +20 @@
    USE phycst         ! NEMO physical constants
    USE sbc_oce
+   USE eosbn2         ! equation of state
 
    USE icb_oce        ! define iceberg arrays
@@ -49 +50 @@
       INTEGER  ::   ii, ij
       REAL(wp) ::   zM, zT, zW, zL, zSST, zVol, zLn, zWn, zTn, znVol, zIC, zDn
+      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
@@ -78 +80 @@
          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%sst, pt%cn, pt%hi, zff, pt%sss )
          !
          zSST = pt%sst
+         zSSS = pt%sss
+         CALL eos_fzp(zSSS,zfzpt)                       ! freezing point
          zIC  = MIN( 1._wp, pt%cn + rn_sicn_shift )     ! Shift sea-ice concentration       !!gm ???
          zM   = pt%mass
@@ -102 +106 @@
 
          ! 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 )
+         zMv = MAX( 7.62d-3*zSST+1.29d-3*(zSST**2)                    , 0._wp ) * z1_rday      ! Buoyant convection at sides (eqn M.A10)
+         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 )
+         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 )
 
          IF( ln_operator_splitting ) THEN      ! Operator split update of volume/mass

NEMO/branches/UKMO/NEMO_4.0.1_ICB_melting_temperature/src/OCE/ICB/icbutl.F90

@@ -74 +74 @@
       ff_e(1:jpi,1:jpj) = ff_f (:,:) 
       tt_e(1:jpi,1:jpj) = sst_m(:,:)
+      ss_e(1:jpi,1:jpj) = sss_m(:,:)
       fr_e(1:jpi,1:jpj) = fr_i (:,:)
       ua_e(1:jpi,1:jpj) = utau (:,:) * umask(:,:,1) ! maybe mask useless because mask applied in sbcblk
@@ -85 +86 @@
       CALL lbc_lnk_icb( 'icbutl', fr_e, 'T', +1._wp, 1, 1 )
       CALL lbc_lnk_icb( 'icbutl', tt_e, 'T', +1._wp, 1, 1 )
+      CALL lbc_lnk_icb( 'icbutl', ss_e, 'T', +1._wp, 1, 1 )
 #if defined key_si3
       hi_e(1:jpi, 1:jpj) = hm_i (:,:)  
@@ -107 +109 @@
    SUBROUTINE icb_utl_interp( pi, pe1, puo, pui, pua, pssh_i,   &
       &                       pj, pe2, pvo, pvi, pva, pssh_j,   &
-      &                       psst, pcn, phi, pff            )
+      &                       psst, pcn, phi, pff, psss        )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE icb_utl_interp  ***
@@ -128 +130 @@
       REAL(wp), INTENT(  out) ::   puo, pvo, pui, pvi, pua, pva   ! ocean, ice and wind speeds
       REAL(wp), INTENT(  out) ::   pssh_i, pssh_j                 ! ssh i- & j-gradients
-      REAL(wp), INTENT(  out) ::   psst, pcn, phi, pff            ! SST, ice concentration, ice thickness, Coriolis
+      REAL(wp), INTENT(  out) ::   psst, pcn, phi, pff, psss      ! SST, ice concentration, ice thickness, Coriolis, SSS
       !
       REAL(wp) ::   zcd, zmod       ! local scalars
@@ -139 +141 @@
       pvo  = icb_utl_bilin_h( vo_e, pi, pj, 'V', .false.  )
       psst = icb_utl_bilin_h( tt_e, pi, pj, 'T', .true.   )    ! SST
+      psss = icb_utl_bilin_h( ss_e, pi, pj, 'T', .true.   )    ! SSS
       pcn  = icb_utl_bilin_h( fr_e, pi, pj, 'T', .true.   )    ! ice concentration
       pff  = icb_utl_bilin_h( ff_e, pi, pj, 'F', .false.  )    ! Coriolis parameter