New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
icbthm.F90 in NEMO/branches/UKMO/NEMO_4.0.4_bouncing_icebergs/src/OCE/ICB – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.4_bouncing_icebergs/src/OCE/ICB/icbthm.F90 @ 14898

Last change on this file since 14898 was 14898, checked in by dancopsey, 16 months ago

Make iceberg top melt a combination of sea ice top melt and ocean solar.

File size: 13.8 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) ::   zqsr, pqsr, zqml, pqml
54      REAL(wp) ::   zSSS, zfzpt
55      REAL(wp) ::   zMv, zMe, zMb, zmelt, zdvo, zdva, zdM, zSs, zdMe, zdMb, zdMv, zMt
56      REAL(wp) ::   zMnew, zMnew1, zMnew2, zheat_hcflux, zheat_latent, z1_12
57      REAL(wp) ::   zMbits, znMbits, zdMbitsE, zdMbitsM, zLbits, zAbits, zMbb
58      REAL(wp) ::   zxi, zyj, zff, z1_rday, z1_e1e2, zdt, z1_dt, z1_dt_e1e2
59      TYPE(iceberg), POINTER ::   this, next
60      TYPE(point)  , POINTER ::   pt
61      !
62      COMPLEX(wp), DIMENSION(jpi,jpj) :: cicb_melt, cicb_hflx
63      !!----------------------------------------------------------------------
64      !
65      !! initialiaze cicb_melt and cicb_heat
66      cicb_melt = CMPLX( 0.e0, 0.e0, wp ) 
67      cicb_hflx = CMPLX( 0.e0, 0.e0, wp ) 
68      !
69      z1_rday = 1._wp / rday
70      z1_12   = 1._wp / 12._wp
71      zdt     = berg_dt
72      z1_dt   = 1._wp / zdt
73      !
74      ! we're either going to ignore berg fresh water melt flux and associated heat
75      ! or we pass it into the ocean, so at this point we set them both to zero,
76      ! accumulate the contributions to them from each iceberg in the while loop following
77      ! and then pass them (or not) to the ocean
78      !
79      berg_grid%floating_melt(:,:) = 0._wp
80      ! calving_hflx re-used here as temporary workspace for the heat flux associated with melting
81      berg_grid%calving_hflx(:,:)  = 0._wp
82      !
83      this => first_berg
84      DO WHILE( ASSOCIATED(this) )
85         !
86         pt => this%current_point
87         nknberg = this%number(1)
88         CALL icb_utl_interp( pt%xi, pt%e1, pt%uo, pt%ui, pt%ua, pt%ssh_x,   &
89            &                 pt%yj, pt%e2, pt%vo, pt%vi, pt%va, pt%ssh_y,   &
90            &                 pt%sst, pt%cn, pt%hi, zff, pt%sss, pqsr, pqml )
91         !
92         zSST = pt%sst
93         zSSS = pt%sss
94         zqsr = MAX( pqsr, 0._wp )                        ! Downward oceanic surface solar flux
95         zqml = MAX( pqml, 0._wp )                        ! Sea ice top melt
96         CALL eos_fzp(zSSS,zfzpt)                       ! freezing point
97         zIC  = MIN( 1._wp, pt%cn + rn_sicn_shift )     ! Shift sea-ice concentration       !!gm ???
98         zM   = pt%mass
99         zT   = pt%thickness                               ! total thickness
100       ! D   = (rn_rho_bergs/pp_rho_seawater)*zT ! draught (keel depth)
101       ! F   = zT - D ! freeboard
102         zW   = pt%width
103         zL   = pt%length
104         zxi  = pt%xi                                      ! position in (i,j) referential
105         zyj  = pt%yj
106         ii  = INT( zxi + 0.5 )                            ! T-cell of the berg
107         ii  = mi1( ii )
108         ij  = INT( zyj + 0.5 )             
109         ij  = mj1( ij )
110         zVol = zT * zW * zL
111
112
113         ! Environment
114         zdvo = SQRT( (pt%uvel-pt%uo)**2 + (pt%vvel-pt%vo)**2 )
115         zdva = SQRT( (pt%ua  -pt%uo)**2 + (pt%va  -pt%vo)**2 )
116         zSs  = 1.5_wp * SQRT( zdva ) + 0.1_wp * zdva                ! Sea state      (eqn M.A9)
117
118         ! Melt rates in m/s (i.e. division by rday)
119         zMv = MAX( 7.62d-3*zSST+1.29d-3*(zSST**2)                    , 0._wp ) * z1_rday      ! Buoyant convection at sides (eqn M.A10)
120         IF ( zSST > zfzpt ) THEN                                                              ! Calculate basal melting only if SST above freezing point 
121            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 )
122         ELSE
123            zMb = 0._wp                                                                        ! No basal melting if SST below freezing point     
124         ENDIF
125         zMe = MAX( z1_12*(zSST+2.)*zSs*(1._wp+COS(rpi*(zIC**3)))     , 0._wp ) * z1_rday      ! Wave erosion                (eqn M.A8 )
126
127         ! Energy involved in top melt is a combination of 10% of solar flux (for portion of grid box that is ocean) and
128         ! 100% of sea ice top melt (for portion of grid box that is sea ice). Apply a maximum amount of melt of 5E-6 m/s
129         ! to keep this term lower than the other terms when the berg is in mid ocean.
130         zMt = MIN( (zqsr*0.1_wp + zqml) / (rLfus * rhoi) , 5.0E-6_wp )
131
132         IF( ln_operator_splitting ) THEN      ! Operator split update of volume/mass
133            zTn    = MAX( zT - zMb*zdt - zMt*zdt , 0._wp )         ! new total thickness (m)
134            znVol  = zTn * zW * zL                       ! new volume (m^3)
135            zMnew1 = ( znVol / zVol ) * zM               ! new mass (kg)
136            zdMb   = zM - zMnew1                         ! mass lost to basal melting (>0) (kg)
137            !
138            zLn    = MAX( zL - zMv*zdt , 0._wp )         ! new length (m)
139            zWn    = MAX( zW - zMv*zdt , 0._wp )         ! new width (m)
140            znVol  = zTn * zWn * zLn                     ! new volume (m^3)
141            zMnew2 = ( znVol / zVol ) * zM               ! new mass (kg)
142            zdMv   = zMnew1 - zMnew2                     ! mass lost to buoyant convection (>0) (kg)
143            !
144            zLn    = MAX( zLn - zMe*zdt , 0._wp )        ! new length (m)
145            zWn    = MAX( zWn - zMe*zdt , 0._wp )        ! new width (m)
146            znVol  = zTn * zWn * zLn                     ! new volume (m^3)
147            zMnew  = ( znVol / zVol ) * zM               ! new mass (kg)
148            zdMe   = zMnew2 - zMnew                      ! mass lost to erosion (>0) (kg)
149            zdM    = zM - zMnew                          ! mass lost to all erosion and melting (>0) (kg)
150            !
151         ELSE                                         ! Update dimensions of berg
152            zLn = MAX( zL -(zMv+zMe)*zdt ,0._wp )        ! (m)
153            zWn = MAX( zW -(zMv+zMe)*zdt ,0._wp )        ! (m)
154            zTn = MAX( zT - zMb*zdt - zMt*zdt ,0._wp )   ! (m)
155            ! Update volume and mass of berg
156            znVol = zTn*zWn*zLn                          ! (m^3)
157            zMnew = (znVol/zVol)*zM                      ! (kg)
158            zdM   = zM - zMnew                           ! (kg)
159            zdMb = (zM/zVol) * (zW*   zL ) *(zMb+zMt)*zdt      ! approx. mass loss to basal melting (kg)
160            zdMe = (zM/zVol) * (zT*(zW+zL)) *zMe*zdt     ! approx. mass lost to erosion (kg)
161            zdMv = (zM/zVol) * (zT*(zW+zL)) *zMv*zdt     ! approx. mass loss to buoyant convection (kg)
162         ENDIF
163
164         IF( rn_bits_erosion_fraction > 0._wp ) THEN     ! Bergy bits
165            !
166            zMbits   = pt%mass_of_bits                                               ! mass of bergy bits (kg)
167            zdMbitsE = rn_bits_erosion_fraction * zdMe                               ! change in mass of bits (kg)
168            znMbits  = zMbits + zdMbitsE                                             ! add new bergy bits to mass (kg)
169            zLbits   = MIN( zL, zW, zT, 40._wp )                                     ! assume bergy bits are smallest dimension or 40 meters
170            zAbits   = ( zMbits / rn_rho_bergs ) / zLbits                            ! Effective bottom area (assuming T=Lbits)
171            zMbb     = MAX( 0.58_wp*(zdvo**0.8_wp)*(zSST+2._wp) /   &
172               &                              ( zLbits**0.2_wp ) , 0._wp ) * z1_rday ! Basal turbulent melting (for bits)
173            zMbb     = rn_rho_bergs * zAbits * zMbb                                  ! in kg/s
174            zdMbitsM = MIN( zMbb*zdt , znMbits )                                     ! bergy bits mass lost to melting (kg)
175            znMbits  = znMbits-zdMbitsM                                              ! remove mass lost to bergy bits melt
176            IF( zMnew <= 0._wp ) THEN                                                ! if parent berg has completely melted then
177               zdMbitsM = zdMbitsM + znMbits                                         ! instantly melt all the bergy bits
178               znMbits  = 0._wp
179            ENDIF
180         ELSE                                                     ! No bergy bits
181            zAbits   = 0._wp
182            zdMbitsE = 0._wp
183            zdMbitsM = 0._wp
184            znMbits  = pt%mass_of_bits                             ! retain previous value incase non-zero
185         ENDIF
186
187
188
189         ! use tmask rather than tmask_i when dealing with icebergs
190         IF( tmask(ii,ij,1) /= 0._wp ) THEN    ! Add melting to the grid and field diagnostics
191            z1_e1e2    = r1_e1e2t(ii,ij) * this%mass_scaling
192            z1_dt_e1e2 = z1_dt * z1_e1e2
193            !
194            ! iceberg melt
195            !! the use of DDPDD function for the cumulative sum is needed for reproducibility
196            zmelt    = ( zdM - ( zdMbitsE - zdMbitsM ) ) * z1_dt   ! kg/s
197            CALL DDPDD( CMPLX( zmelt * z1_e1e2, 0.e0, wp ), cicb_melt(ii,ij) )
198            !
199            ! iceberg heat flux
200            !! the use of DDPDD function for the cumulative sum is needed for reproducibility
201            !! NB. The src_calving_hflx field is currently hardwired to zero in icb_stp, which means that the
202            !!     heat density of the icebergs is zero and the heat content flux to the ocean from iceberg
203            !!     melting is always zero. Leaving the term in the code until such a time as this is fixed. DS.
204            zheat_hcflux = zmelt * pt%heat_density       ! heat content flux : kg/s x J/kg = J/s
205            zheat_latent = - zmelt * rLfus               ! latent heat flux:  kg/s x J/kg = J/s
206            CALL DDPDD( CMPLX( ( zheat_hcflux + zheat_latent ) * z1_e1e2, 0.e0, wp ), cicb_hflx(ii,ij) )
207            !
208            ! diagnostics
209            CALL icb_dia_melt( ii, ij, zMnew, zheat_hcflux, zheat_latent, this%mass_scaling,       &
210               &                       zdM, zdMbitsE, zdMbitsM, zdMb, zdMe,   &
211               &                       zdMv, z1_dt_e1e2 )
212         ELSE
213            WRITE(numout,*) 'icb_thm: berg ',this%number(:),' appears to have grounded  at ',narea,ii,ij
214            CALL icb_utl_print_berg( this, kt )
215            WRITE(numout,*) 'msk=',tmask(ii,ij,1), e1e2t(ii,ij)
216            CALL ctl_stop('icb_thm', 'berg appears to have grounded!')
217         ENDIF
218
219         ! Rolling
220         zDn = ( rn_rho_bergs / pp_rho_seawater ) * zTn       ! draught (keel depth)
221         IF( zDn > 0._wp .AND. MAX(zWn,zLn) < SQRT( 0.92*(zDn**2) + 58.32*zDn ) ) THEN
222            zT  = zTn
223            zTn = zWn
224            zWn = zT
225         ENDIF
226
227         ! Store the new state of iceberg (with L>W)
228         pt%mass         = zMnew
229         pt%mass_of_bits = znMbits
230         pt%thickness    = zTn
231         pt%width        = MIN( zWn , zLn )
232         pt%length       = MAX( zWn , zLn )
233
234         next=>this%next
235
236!!gm  add a test to avoid over melting ?
237
238         IF( zMnew <= 0._wp ) THEN       ! Delete the berg if completely melted
239            CALL icb_utl_delete( first_berg, this )
240            !
241         ELSE                            ! Diagnose mass distribution on grid
242            z1_e1e2 = r1_e1e2t(ii,ij) * this%mass_scaling
243            CALL icb_dia_size( ii, ij, zWn, zLn, zAbits,   &
244               &               this%mass_scaling, zMnew, znMbits, z1_e1e2 )
245         ENDIF
246         !
247         this=>next
248         !
249      END DO
250      !
251      berg_grid%floating_melt = REAL(cicb_melt,wp)    ! kg/m2/s
252      berg_grid%calving_hflx  = REAL(cicb_hflx,wp)
253      !
254      ! now use melt and associated heat flux in ocean (or not)
255      !
256      IF(.NOT. ln_passive_mode ) THEN
257         emp (:,:) = emp (:,:) - berg_grid%floating_melt(:,:)
258         qns (:,:) = qns (:,:) + berg_grid%calving_hflx (:,:) 
259      ENDIF
260      !
261   END SUBROUTINE icb_thm
262
263   !!======================================================================
264END MODULE icbthm
Note: See TracBrowser for help on using the repository browser.