source: NEMO/trunk/src/OCE/ICB/icbthm.F90 @ 13281

Last change on this file since 13281 was 13281, checked in by ayoung, 2 months ago

Merging changes for ticket #2407. Patch for iceberg melting stability. Fully SETTE tested: this merge introduces a change to the results in ORCA2_ICE_PISCES.

  • Property svn:keywords set to Id
File size: 13.2 KB
Line 
1MODULE icbthm
2   !!======================================================================
3   !!                       ***  MODULE  icbthm  ***
4   !! Icebergs:  thermodynamics routines for icebergs
5   !!======================================================================
6   !! History : 3.3.1 !  2010-01  (Martin&Adcroft) Original code
7   !!            -    !  2011-03  (Madec)          Part conversion to NEMO form
8   !!            -    !                            Removal of mapping from another grid
9   !!            -    !  2011-04  (Alderson)       Split into separate modules
10   !!            -    !  2011-05  (Alderson)       Use tmask instead of tmask_i
11   !!----------------------------------------------------------------------
12   !!----------------------------------------------------------------------
13   !!   icb_thm : initialise
14   !!             reference for equations - M = Martin + Adcroft, OM 34, 2010
15   !!----------------------------------------------------------------------
16   USE par_oce        ! NEMO parameters
17   USE dom_oce        ! NEMO domain
18   USE in_out_manager ! NEMO IO routines, numout in particular
19   USE lib_mpp        ! NEMO MPI routines, ctl_stop in particular
20   USE phycst         ! NEMO physical constants
21   USE sbc_oce
22   USE eosbn2         ! equation of state
23   USE lib_fortran, ONLY : DDPDD
24
25   USE icb_oce        ! define iceberg arrays
26   USE icbutl         ! iceberg utility routines
27   USE icbdia         ! iceberg budget routines
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   icb_thm ! routine called in icbstp.F90 module
33
34   !!----------------------------------------------------------------------
35   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
36   !! $Id$
37   !! Software governed by the CeCILL license (see ./LICENSE)
38   !!----------------------------------------------------------------------
39CONTAINS
40
41   SUBROUTINE icb_thm( kt )
42      !!----------------------------------------------------------------------
43      !!                  ***  ROUTINE icb_thm  ***
44      !!
45      !! ** Purpose :   compute the iceberg thermodynamics.
46      !!
47      !! ** Method  : - See Martin & Adcroft, Ocean Modelling 34, 2010
48      !!----------------------------------------------------------------------
49      INTEGER, INTENT(in) ::   kt   ! timestep number, just passed to icb_utl_print_berg
50      !
51      INTEGER  ::   ii, ij
52      REAL(wp) ::   zM, zT, zW, zL, zSST, zVol, zLn, zWn, zTn, znVol, zIC, zDn
53      REAL(wp) ::   zSSS, zfzpt
54      REAL(wp) ::   zMv, zMe, zMb, zmelt, zdvo, zdva, zdM, zSs, zdMe, zdMb, zdMv
55      REAL(wp) ::   zMnew, zMnew1, zMnew2, zheat_hcflux, zheat_latent, z1_12
56      REAL(wp) ::   zMbits, znMbits, zdMbitsE, zdMbitsM, zLbits, zAbits, zMbb
57      REAL(wp) ::   zxi, zyj, zff, z1_rday, z1_e1e2, zdt, z1_dt, z1_dt_e1e2
58      TYPE(iceberg), POINTER ::   this, next
59      TYPE(point)  , POINTER ::   pt
60      !
61      COMPLEX(dp), DIMENSION(jpi,jpj) :: cicb_melt, cicb_hflx
62      !!----------------------------------------------------------------------
63      !
64      !! initialiaze cicb_melt and cicb_heat
65      cicb_melt = CMPLX( 0.e0, 0.e0, dp ) 
66      cicb_hflx = CMPLX( 0.e0, 0.e0, dp ) 
67      !
68      z1_rday = 1._wp / rday
69      z1_12   = 1._wp / 12._wp
70      zdt     = berg_dt
71      z1_dt   = 1._wp / zdt
72      !
73      ! we're either going to ignore berg fresh water melt flux and associated heat
74      ! or we pass it into the ocean, so at this point we set them both to zero,
75      ! accumulate the contributions to them from each iceberg in the while loop following
76      ! and then pass them (or not) to the ocean
77      !
78      berg_grid%floating_melt(:,:) = 0._wp
79      ! calving_hflx re-used here as temporary workspace for the heat flux associated with melting
80      berg_grid%calving_hflx(:,:)  = 0._wp
81      !
82      this => first_berg
83      DO WHILE( ASSOCIATED(this) )
84         !
85         pt => this%current_point
86         nknberg = this%number(1)
87         CALL icb_utl_interp( pt%xi, pt%e1, pt%uo, pt%ui, pt%ua, pt%ssh_x,   &
88            &                 pt%yj, pt%e2, pt%vo, pt%vi, pt%va, pt%ssh_y,   &
89            &                 pt%sst, pt%cn, pt%hi, zff, pt%sss )
90         !
91         zSST = pt%sst
92         zSSS = pt%sss
93         CALL eos_fzp(zSSS,zfzpt)                       ! freezing point
94         zIC  = MIN( 1._wp, pt%cn + rn_sicn_shift )     ! Shift sea-ice concentration       !!gm ???
95         zM   = pt%mass
96         zT   = pt%thickness                               ! total thickness
97       ! D   = (rn_rho_bergs/pp_rho_seawater)*zT ! draught (keel depth)
98       ! F   = zT - D ! freeboard
99         zW   = pt%width
100         zL   = pt%length
101         zxi  = pt%xi                                      ! position in (i,j) referential
102         zyj  = pt%yj
103         ii  = INT( zxi + 0.5 )                            ! T-cell of the berg
104         ii  = mi1( ii )
105         ij  = INT( zyj + 0.5 )             
106         ij  = mj1( ij )
107         zVol = zT * zW * zL
108
109         ! Environment
110         zdvo = SQRT( (pt%uvel-pt%uo)**2 + (pt%vvel-pt%vo)**2 )
111         zdva = SQRT( (pt%ua  -pt%uo)**2 + (pt%va  -pt%vo)**2 )
112         zSs  = 1.5_wp * SQRT( zdva ) + 0.1_wp * zdva                ! Sea state      (eqn M.A9)
113
114         ! Melt rates in m/s (i.e. division by rday)
115         zMv = MAX( 7.62d-3*zSST+1.29d-3*(zSST**2)                    , 0._wp ) * z1_rday      ! Buoyant convection at sides (eqn M.A10)
116         IF ( zSST > zfzpt ) THEN                                                              ! Calculate basal melting only if SST above freezing point 
117            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 )
118         ELSE
119            zMb = 0._wp                                                                        ! No basal melting if SST below freezing point     
120         ENDIF
121         zMe = MAX( z1_12*(zSST+2.)*zSs*(1._wp+COS(rpi*(zIC**3)))     , 0._wp ) * z1_rday      ! Wave erosion                (eqn M.A8 )
122
123         IF( ln_operator_splitting ) THEN      ! Operator split update of volume/mass
124            zTn    = MAX( zT - zMb*zdt , 0._wp )         ! new total thickness (m)
125            znVol  = zTn * zW * zL                       ! new volume (m^3)
126            zMnew1 = ( znVol / zVol ) * zM               ! new mass (kg)
127            zdMb   = zM - zMnew1                         ! mass lost to basal melting (>0) (kg)
128            !
129            zLn    = MAX( zL - zMv*zdt , 0._wp )         ! new length (m)
130            zWn    = MAX( zW - zMv*zdt , 0._wp )         ! new width (m)
131            znVol  = zTn * zWn * zLn                     ! new volume (m^3)
132            zMnew2 = ( znVol / zVol ) * zM               ! new mass (kg)
133            zdMv   = zMnew1 - zMnew2                     ! mass lost to buoyant convection (>0) (kg)
134            !
135            zLn    = MAX( zLn - zMe*zdt , 0._wp )        ! new length (m)
136            zWn    = MAX( zWn - zMe*zdt , 0._wp )        ! new width (m)
137            znVol  = zTn * zWn * zLn                     ! new volume (m^3)
138            zMnew  = ( znVol / zVol ) * zM               ! new mass (kg)
139            zdMe   = zMnew2 - zMnew                      ! mass lost to erosion (>0) (kg)
140            zdM    = zM - zMnew                          ! mass lost to all erosion and melting (>0) (kg)
141            !
142         ELSE                                         ! Update dimensions of berg
143            zLn = MAX( zL -(zMv+zMe)*zdt ,0._wp )        ! (m)
144            zWn = MAX( zW -(zMv+zMe)*zdt ,0._wp )        ! (m)
145            zTn = MAX( zT - zMb    *zdt ,0._wp )         ! (m)
146            ! Update volume and mass of berg
147            znVol = zTn*zWn*zLn                          ! (m^3)
148            zMnew = (znVol/zVol)*zM                      ! (kg)
149            zdM   = zM - zMnew                           ! (kg)
150            zdMb = (zM/zVol) * (zW*   zL ) *zMb*zdt      ! approx. mass loss to basal melting (kg)
151            zdMe = (zM/zVol) * (zT*(zW+zL)) *zMe*zdt     ! approx. mass lost to erosion (kg)
152            zdMv = (zM/zVol) * (zT*(zW+zL)) *zMv*zdt     ! approx. mass loss to buoyant convection (kg)
153         ENDIF
154
155         IF( rn_bits_erosion_fraction > 0._wp ) THEN     ! Bergy bits
156            !
157            zMbits   = pt%mass_of_bits                                               ! mass of bergy bits (kg)
158            zdMbitsE = rn_bits_erosion_fraction * zdMe                               ! change in mass of bits (kg)
159            znMbits  = zMbits + zdMbitsE                                             ! add new bergy bits to mass (kg)
160            zLbits   = MIN( zL, zW, zT, 40._wp )                                     ! assume bergy bits are smallest dimension or 40 meters
161            zAbits   = ( zMbits / rn_rho_bergs ) / zLbits                            ! Effective bottom area (assuming T=Lbits)
162            zMbb     = MAX( 0.58_wp*(zdvo**0.8_wp)*(zSST+2._wp) /   &
163               &                              ( zLbits**0.2_wp ) , 0._wp ) * z1_rday ! Basal turbulent melting (for bits)
164            zMbb     = rn_rho_bergs * zAbits * zMbb                                  ! in kg/s
165            zdMbitsM = MIN( zMbb*zdt , znMbits )                                     ! bergy bits mass lost to melting (kg)
166            znMbits  = znMbits-zdMbitsM                                              ! remove mass lost to bergy bits melt
167            IF( zMnew == 0._wp ) THEN                                                ! if parent berg has completely melted then
168               zdMbitsM = zdMbitsM + znMbits                                         ! instantly melt all the bergy bits
169               znMbits  = 0._wp
170            ENDIF
171         ELSE                                                     ! No bergy bits
172            zAbits   = 0._wp
173            zdMbitsE = 0._wp
174            zdMbitsM = 0._wp
175            znMbits  = pt%mass_of_bits                             ! retain previous value incase non-zero
176         ENDIF
177
178         ! use tmask rather than tmask_i when dealing with icebergs
179         IF( tmask(ii,ij,1) /= 0._wp ) THEN    ! Add melting to the grid and field diagnostics
180            z1_e1e2    = r1_e1e2t(ii,ij) * this%mass_scaling
181            z1_dt_e1e2 = z1_dt * z1_e1e2
182            !
183            ! iceberg melt
184            !! the use of DDPDD function for the cumulative sum is needed for reproducibility
185            zmelt    = ( zdM - ( zdMbitsE - zdMbitsM ) ) * z1_dt   ! kg/s
186            CALL DDPDD( CMPLX( zmelt * z1_e1e2, 0.e0, dp ), cicb_melt(ii,ij) )
187            !
188            ! iceberg heat flux
189            !! the use of DDPDD function for the cumulative sum is needed for reproducibility
190            !! NB. The src_calving_hflx field is currently hardwired to zero in icb_stp, which means that the
191            !!     heat density of the icebergs is zero and the heat content flux to the ocean from iceberg
192            !!     melting is always zero. Leaving the term in the code until such a time as this is fixed. DS.
193            zheat_hcflux = zmelt * pt%heat_density       ! heat content flux : kg/s x J/kg = J/s
194            zheat_latent = - zmelt * rLfus               ! latent heat flux:  kg/s x J/kg = J/s
195            CALL DDPDD( CMPLX( ( zheat_hcflux + zheat_latent ) * z1_e1e2, 0.e0, dp ), cicb_hflx(ii,ij) )
196            !
197            ! diagnostics
198            CALL icb_dia_melt( ii, ij, zMnew, zheat_hcflux, zheat_latent, this%mass_scaling,       &
199               &                       zdM, zdMbitsE, zdMbitsM, zdMb, zdMe,   &
200               &                       zdMv, z1_dt_e1e2 )
201         ELSE
202            WRITE(numout,*) 'icb_thm: berg ',this%number(:),' appears to have grounded  at ',narea,ii,ij
203            CALL icb_utl_print_berg( this, kt )
204            WRITE(numout,*) 'msk=',tmask(ii,ij,1), e1e2t(ii,ij)
205            CALL ctl_stop('icb_thm', 'berg appears to have grounded!')
206         ENDIF
207
208         ! Rolling
209         zDn = ( rn_rho_bergs / pp_rho_seawater ) * zTn       ! draught (keel depth)
210         IF( zDn > 0._wp .AND. MAX(zWn,zLn) < SQRT( 0.92*(zDn**2) + 58.32*zDn ) ) THEN
211            zT  = zTn
212            zTn = zWn
213            zWn = zT
214         ENDIF
215
216         ! Store the new state of iceberg (with L>W)
217         pt%mass         = zMnew
218         pt%mass_of_bits = znMbits
219         pt%thickness    = zTn
220         pt%width        = MIN( zWn , zLn )
221         pt%length       = MAX( zWn , zLn )
222
223         next=>this%next
224
225!!gm  add a test to avoid over melting ?
226
227         IF( zMnew <= 0._wp ) THEN       ! Delete the berg if completely melted
228            CALL icb_utl_delete( first_berg, this )
229            !
230         ELSE                            ! Diagnose mass distribution on grid
231            z1_e1e2 = r1_e1e2t(ii,ij) * this%mass_scaling
232            CALL icb_dia_size( ii, ij, zWn, zLn, zAbits,   &
233               &               this%mass_scaling, zMnew, znMbits, z1_e1e2 )
234         ENDIF
235         !
236         this=>next
237         !
238      END DO
239      !
240      berg_grid%floating_melt = REAL(cicb_melt,dp)    ! kg/m2/s
241      berg_grid%calving_hflx  = REAL(cicb_hflx,dp)
242      !
243      ! now use melt and associated heat flux in ocean (or not)
244      !
245      IF(.NOT. ln_passive_mode ) THEN
246         emp (:,:) = emp (:,:) - berg_grid%floating_melt(:,:)
247         qns (:,:) = qns (:,:) + berg_grid%calving_hflx (:,:) 
248      ENDIF
249      !
250   END SUBROUTINE icb_thm
251
252   !!======================================================================
253END MODULE icbthm
Note: See TracBrowser for help on using the repository browser.