source: NEMO/branches/UKMO/NEMO_4.0.1_ICB_melting_temperature/src/OCE/ICB/icbthm.F90 @ 11897

Last change on this file since 11897 was 11897, checked in by cguiavarch, 12 months ago

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

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