source: NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbthm.F90 @ 13357

Last change on this file since 13357 was 13357, checked in by mathiot, 8 weeks ago

ticket #1900: add changes for Merino 2016 works. Results unchanged if flag ln_M2016 is set to .false. . Remaining step: test ln_M2016 = .true.

  • Property svn:keywords set to Id
File size: 14.4 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 lib_fortran, ONLY : DDPDD
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, jk, ikb
51      REAL(wp) ::   zM, zT, zW, zL, zSST, zVol, zLn, zWn, zTn, znVol, zIC, zDn, zD, zvb, zub, ztb
52      REAL(wp) ::   zMv, zMe, zMb, zmelt, zdvo, zdvob, zdva, zdM, zSs, zdMe, zdMb, zdMv
53      REAL(wp) ::   zMnew, zMnew1, zMnew2, zheat_hcflux, zheat_latent, z1_12
54      REAL(wp) ::   zMbits, znMbits, zdMbitsE, zdMbitsM, zLbits, zAbits, zMbb
55      REAL(wp) ::   zxi, zyj, zff, z1_rday, z1_e1e2, zdt, z1_dt, z1_dt_e1e2, zdepw
56      REAL(wp), DIMENSION(jpk) :: ztoce, zuoce, zvoce, ze3t, ztmp
57      TYPE(iceberg), POINTER ::   this, next
58      TYPE(point)  , POINTER ::   pt
59      !
60      COMPLEX(dp), DIMENSION(jpi,jpj) :: cicb_melt, cicb_hflx
61      !!----------------------------------------------------------------------
62      !
63      !! initialiaze cicb_melt and cicb_heat
64      cicb_melt = CMPLX( 0.e0, 0.e0, dp ) 
65      cicb_hflx = CMPLX( 0.e0, 0.e0, dp ) 
66      !
67      z1_rday = 1._wp / rday
68      z1_12   = 1._wp / 12._wp
69      zdt     = berg_dt
70      z1_dt   = 1._wp / zdt
71      !
72      ! we're either going to ignore berg fresh water melt flux and associated heat
73      ! or we pass it into the ocean, so at this point we set them both to zero,
74      ! accumulate the contributions to them from each iceberg in the while loop following
75      ! and then pass them (or not) to the ocean
76      !
77      berg_grid%floating_melt(:,:) = 0._wp
78      ! calving_hflx re-used here as temporary workspace for the heat flux associated with melting
79      berg_grid%calving_hflx(:,:)  = 0._wp
80      !
81      this => first_berg
82      DO WHILE( ASSOCIATED(this) )
83         !
84         pt => this%current_point
85         nknberg = this%number(1)
86
87         CALL icb_utl_interp( pt%xi, pt%yj,           &   ! position
88             &                 pssu=pt%ssu, pua=pt%ua, &   ! oce/atm velocities
89             &                 pssv=pt%ssv, pva=pt%va, &   ! oce/atm velocities
90             &                 psst=pt%sst, pcn=pt%cn )
91
92         IF ( nn_sample_rate > 0 .AND. MOD(kt-1,nn_sample_rate) == 0 ) THEN
93            CALL icb_utl_interp( pt%xi, pt%yj, pe1=pt%e1, pe2=pt%e2,                 &
94               &                 pui=pt%ui, pssh_i=pt%ssh_x, &
95               &                 pvi=pt%vi, pssh_j=pt%ssh_y, &
96               &                 phi=pt%hi,                  &
97               &                 plat=pt%lat, plon=pt%lon )
98         END IF
99         !
100         zSST = pt%sst
101         zIC  = MIN( 1._wp, pt%cn + rn_sicn_shift )     ! Shift sea-ice concentration       !!gm ???
102         zM   = pt%mass
103         zT   = pt%thickness                               ! total thickness
104         zD   = rho_berg_1_oce * zT          ! draught (keel depth)
105       ! F   = zT - D ! freeboard
106         zW   = pt%width
107         zL   = pt%length
108         zxi  = pt%xi                                      ! position in (i,j) referential
109         zyj  = pt%yj
110         ii  = INT( zxi + 0.5 )                            ! T-cell of the berg
111         ii  = mi1( ii )
112         ij  = INT( zyj + 0.5 )             
113         ij  = mj1( ij )
114         zVol = zT * zW * zL
115
116         ! Environment
117         ! default sst, ssu and ssv
118         ! ln_M2016: use temp, u and v profile
119         IF ( ln_M2016 ) THEN
120
121            ! load t, u, v and e3 profile at icb position
122            CALL icb_utl_interp( pt%xi, pt%yj, ptoce=ztoce, puoce=zuoce, pvoce=zvoce, pe3t=ze3t )
123           
124            !compute bottom level
125            CALL icb_utl_getkb( pt%kb, ze3t, zD )
126
127            ikb = MIN(pt%kb,mbkt(ii,ij))
128            ztb = ztoce(ikb)                                                ! basal temperature
129            zub = zuoce(ikb)
130            zvb = zvoce(ikb)
131         ELSE
132            ztb = pt%sst
133            zub = pt%ssu
134            zvb = pt%ssv
135         END IF
136
137         zdvob = SQRT( (pt%uvel-zub)**2 + (pt%vvel-zvb)**2 )        ! relative basal velocity
138         zdva  = SQRT( (pt%ua  -pt%ssu)**2 + (pt%va  -pt%ssv)**2 )  ! relative wind
139         zSs   = 1.5_wp * SQRT( zdva ) + 0.1_wp * zdva              ! Sea state      (eqn M.A9)
140         !
141         ! Melt rates in m/s (i.e. division by rday)
142         IF ( ln_M2016 ) THEN
143            ! Buoyant convection at sides (eqn M.A10) but averaging along all the iceberg draft
144            ztmp(:) = MAX( 7.62d-3*ztoce(:)+1.29d-3*(ztoce(:)**2), 0._wp ) * z1_rday
145            zMv = 0.0; zdepw = 0.0
146            DO jk = 1,ikb-1
147               zMv = zMv + ze3t(jk)*ztmp(jk)
148               zdepw = zdepw + ze3t(jk)
149            END DO
150            zMv = (zMv + (zD - zdepw)*ztmp(ikb)) / zD
151         ELSE
152            zMv = MAX( 7.62d-3*zSST+1.29d-3*(zSST**2), 0._wp ) * z1_rday
153         END IF
154
155         ! Basal turbulent melting     (eqn M.A7 ) but using the basal temperature
156         zMb = MAX( 0.58_wp*(zdvob**0.8_wp)*(ztb+4.0_wp)/(zL**0.2_wp) , 0._wp ) * z1_rday
157
158         ! Wave erosion                (eqn M.A8 )
159         zMe = MAX( z1_12*(zSST+2.)*zSs*(1._wp+COS(rpi*(zIC**3)))     , 0._wp ) * z1_rday
160
161         IF( ln_operator_splitting ) THEN      ! Operator split update of volume/mass
162            zTn    = MAX( zT - zMb*zdt , 0._wp )         ! new total thickness (m)
163            znVol  = zTn * zW * zL                       ! new volume (m^3)
164            zMnew1 = ( znVol / zVol ) * zM               ! new mass (kg)
165            zdMb   = zM - zMnew1                         ! mass lost to basal melting (>0) (kg)
166            !
167            zLn    = MAX( zL - zMv*zdt , 0._wp )         ! new length (m)
168            zWn    = MAX( zW - zMv*zdt , 0._wp )         ! new width (m)
169            znVol  = zTn * zWn * zLn                     ! new volume (m^3)
170            zMnew2 = ( znVol / zVol ) * zM               ! new mass (kg)
171            zdMv   = zMnew1 - zMnew2                     ! mass lost to buoyant convection (>0) (kg)
172            !
173            zLn    = MAX( zLn - zMe*zdt , 0._wp )        ! new length (m)
174            zWn    = MAX( zWn - zMe*zdt , 0._wp )        ! new width (m)
175            znVol  = zTn * zWn * zLn                     ! new volume (m^3)
176            zMnew  = ( znVol / zVol ) * zM               ! new mass (kg)
177            zdMe   = zMnew2 - zMnew                      ! mass lost to erosion (>0) (kg)
178            zdM    = zM - zMnew                          ! mass lost to all erosion and melting (>0) (kg)
179            !
180         ELSE                                         ! Update dimensions of berg
181            zLn = MAX( zL -(zMv+zMe)*zdt ,0._wp )        ! (m)
182            zWn = MAX( zW -(zMv+zMe)*zdt ,0._wp )        ! (m)
183            zTn = MAX( zT - zMb    *zdt ,0._wp )         ! (m)
184            ! Update volume and mass of berg
185            znVol = zTn*zWn*zLn                          ! (m^3)
186            zMnew = (znVol/zVol)*zM                      ! (kg)
187            zdM   = zM - zMnew                           ! (kg)
188            zdMb = (zM/zVol) * (zW*   zL ) *zMb*zdt      ! approx. mass loss to basal melting (kg)
189            zdMe = (zM/zVol) * (zT*(zW+zL)) *zMe*zdt     ! approx. mass lost to erosion (kg)
190            zdMv = (zM/zVol) * (zT*(zW+zL)) *zMv*zdt     ! approx. mass loss to buoyant convection (kg)
191         ENDIF
192
193         IF( rn_bits_erosion_fraction > 0._wp ) THEN     ! Bergy bits
194            !
195            zMbits   = pt%mass_of_bits                                               ! mass of bergy bits (kg)
196            zdMbitsE = rn_bits_erosion_fraction * zdMe                               ! change in mass of bits (kg)
197            znMbits  = zMbits + zdMbitsE                                             ! add new bergy bits to mass (kg)
198            zLbits   = MIN( zL, zW, zT, 40._wp )                                     ! assume bergy bits are smallest dimension or 40 meters
199            zAbits   = ( zMbits / rn_rho_bergs ) / zLbits                            ! Effective bottom area (assuming T=Lbits)
200            zMbb     = MAX( 0.58_wp*(zdvo**0.8_wp)*(zSST+2._wp) /   &
201               &                              ( zLbits**0.2_wp ) , 0._wp ) * z1_rday ! Basal turbulent melting (for bits)
202            zMbb     = rn_rho_bergs * zAbits * zMbb                                  ! in kg/s
203            zdMbitsM = MIN( zMbb*zdt , znMbits )                                     ! bergy bits mass lost to melting (kg)
204            znMbits  = znMbits-zdMbitsM                                              ! remove mass lost to bergy bits melt
205            IF( zMnew == 0._wp ) THEN                                                ! if parent berg has completely melted then
206               zdMbitsM = zdMbitsM + znMbits                                         ! instantly melt all the bergy bits
207               znMbits  = 0._wp
208            ENDIF
209         ELSE                                                     ! No bergy bits
210            zAbits   = 0._wp
211            zdMbitsE = 0._wp
212            zdMbitsM = 0._wp
213            znMbits  = pt%mass_of_bits                             ! retain previous value incase non-zero
214         ENDIF
215
216         ! use tmask rather than tmask_i when dealing with icebergs
217         IF( tmask(ii,ij,1) /= 0._wp ) THEN    ! Add melting to the grid and field diagnostics
218            z1_e1e2    = r1_e1e2t(ii,ij) * this%mass_scaling
219            z1_dt_e1e2 = z1_dt * z1_e1e2
220            !
221            ! iceberg melt
222            !! the use of DDPDD function for the cumulative sum is needed for reproducibility
223            zmelt    = ( zdM - ( zdMbitsE - zdMbitsM ) ) * z1_dt   ! kg/s
224            CALL DDPDD( CMPLX( zmelt * z1_e1e2, 0.e0, dp ), cicb_melt(ii,ij) )
225            !
226            ! iceberg heat flux
227            !! the use of DDPDD function for the cumulative sum is needed for reproducibility
228            !! NB. The src_calving_hflx field is currently hardwired to zero in icb_stp, which means that the
229            !!     heat density of the icebergs is zero and the heat content flux to the ocean from iceberg
230            !!     melting is always zero. Leaving the term in the code until such a time as this is fixed. DS.
231            zheat_hcflux = zmelt * pt%heat_density       ! heat content flux : kg/s x J/kg = J/s
232            zheat_latent = - zmelt * rLfus               ! latent heat flux:  kg/s x J/kg = J/s
233            CALL DDPDD( CMPLX( ( zheat_hcflux + zheat_latent ) * z1_e1e2, 0.e0, dp ), cicb_hflx(ii,ij) )
234            !
235            ! diagnostics
236            CALL icb_dia_melt( ii, ij, zMnew, zheat_hcflux, zheat_latent, this%mass_scaling,       &
237               &                       zdM, zdMbitsE, zdMbitsM, zdMb, zdMe,   &
238               &                       zdMv, z1_dt_e1e2 )
239         ELSE
240            WRITE(numout,*) 'icb_thm: berg ',this%number(:),' appears to have grounded  at ',narea,ii,ij
241            CALL icb_utl_print_berg( this, kt )
242            WRITE(numout,*) 'msk=',tmask(ii,ij,1), e1e2t(ii,ij)
243            CALL ctl_stop('icb_thm', 'berg appears to have grounded!')
244         ENDIF
245
246         ! Rolling
247         zDn = rho_berg_1_oce * zTn       ! draught (keel depth)
248         IF( zDn > 0._wp .AND. MAX(zWn,zLn) < SQRT( 0.92*(zDn**2) + 58.32*zDn ) ) THEN
249            zT  = zTn
250            zTn = zWn
251            zWn = zT
252         ENDIF
253
254         ! Store the new state of iceberg (with L>W)
255         pt%mass         = zMnew
256         pt%mass_of_bits = znMbits
257         pt%thickness    = zTn
258         pt%width        = MIN( zWn , zLn )
259         pt%length       = MAX( zWn , zLn )
260
261         next=>this%next
262
263!!gm  add a test to avoid over melting ?
264
265         IF( zMnew <= 0._wp ) THEN       ! Delete the berg if completely melted
266            CALL icb_utl_delete( first_berg, this )
267            !
268         ELSE                            ! Diagnose mass distribution on grid
269            z1_e1e2 = r1_e1e2t(ii,ij) * this%mass_scaling
270            CALL icb_dia_size( ii, ij, zWn, zLn, zAbits,   &
271               &               this%mass_scaling, zMnew, znMbits, z1_e1e2 )
272         ENDIF
273         !
274         this=>next
275         !
276      END DO
277      !
278      berg_grid%floating_melt = REAL(cicb_melt,dp)    ! kg/m2/s
279      berg_grid%calving_hflx  = REAL(cicb_hflx,dp)
280      !
281      ! now use melt and associated heat flux in ocean (or not)
282      !
283      IF(.NOT. ln_passive_mode ) THEN
284         emp (:,:) = emp (:,:) - berg_grid%floating_melt(:,:)
285         qns (:,:) = qns (:,:) + berg_grid%calving_hflx (:,:) 
286      ENDIF
287      !
288   END SUBROUTINE icb_thm
289
290   !!======================================================================
291END MODULE icbthm
Note: See TracBrowser for help on using the repository browser.