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

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

ticket #1900: add subroutine icb_utl_zavg, fix issue in icb_utl_getkb, add subroutine test_icb_utl_getkb to test awkward cases (will be removed after the review), simplify icbdyn.F90 icbthm.F90 according to the new routine created.

  • Property svn:keywords set to Id
File size: 14.3 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, zzMv
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            zzMv(:) = MAX( 7.62d-3*ztoce(:)+1.29d-3*(ztoce(:)**2), 0._wp ) * z1_rday
145            CALL icb_utl_zavg(zMv, zzMv, ze3t, zD, ikb )
146         ELSE
147            zMv = MAX( 7.62d-3*zSST+1.29d-3*(zSST**2), 0._wp ) * z1_rday
148         END IF
149
150         ! Basal turbulent melting     (eqn M.A7 ) but using the basal temperature
151         zMb = MAX( 0.58_wp*(zdvob**0.8_wp)*(ztb+4.0_wp)/(zL**0.2_wp) , 0._wp ) * z1_rday
152
153         ! Wave erosion                (eqn M.A8 )
154         zMe = MAX( z1_12*(zSST+2.)*zSs*(1._wp+COS(rpi*(zIC**3)))     , 0._wp ) * z1_rday
155
156         IF( ln_operator_splitting ) THEN      ! Operator split update of volume/mass
157            zTn    = MAX( zT - zMb*zdt , 0._wp )         ! new total thickness (m)
158            znVol  = zTn * zW * zL                       ! new volume (m^3)
159            zMnew1 = ( znVol / zVol ) * zM               ! new mass (kg)
160            zdMb   = zM - zMnew1                         ! mass lost to basal melting (>0) (kg)
161            !
162            zLn    = MAX( zL - zMv*zdt , 0._wp )         ! new length (m)
163            zWn    = MAX( zW - zMv*zdt , 0._wp )         ! new width (m)
164            znVol  = zTn * zWn * zLn                     ! new volume (m^3)
165            zMnew2 = ( znVol / zVol ) * zM               ! new mass (kg)
166            zdMv   = zMnew1 - zMnew2                     ! mass lost to buoyant convection (>0) (kg)
167            !
168            zLn    = MAX( zLn - zMe*zdt , 0._wp )        ! new length (m)
169            zWn    = MAX( zWn - zMe*zdt , 0._wp )        ! new width (m)
170            znVol  = zTn * zWn * zLn                     ! new volume (m^3)
171            zMnew  = ( znVol / zVol ) * zM               ! new mass (kg)
172            zdMe   = zMnew2 - zMnew                      ! mass lost to erosion (>0) (kg)
173            zdM    = zM - zMnew                          ! mass lost to all erosion and melting (>0) (kg)
174            !
175         ELSE                                         ! Update dimensions of berg
176            zLn = MAX( zL -(zMv+zMe)*zdt ,0._wp )        ! (m)
177            zWn = MAX( zW -(zMv+zMe)*zdt ,0._wp )        ! (m)
178            zTn = MAX( zT - zMb    *zdt ,0._wp )         ! (m)
179            ! Update volume and mass of berg
180            znVol = zTn*zWn*zLn                          ! (m^3)
181            zMnew = (znVol/zVol)*zM                      ! (kg)
182            zdM   = zM - zMnew                           ! (kg)
183            zdMb = (zM/zVol) * (zW*   zL ) *zMb*zdt      ! approx. mass loss to basal melting (kg)
184            zdMe = (zM/zVol) * (zT*(zW+zL)) *zMe*zdt     ! approx. mass lost to erosion (kg)
185            zdMv = (zM/zVol) * (zT*(zW+zL)) *zMv*zdt     ! approx. mass loss to buoyant convection (kg)
186         ENDIF
187
188         IF( rn_bits_erosion_fraction > 0._wp ) THEN     ! Bergy bits
189            !
190            zMbits   = pt%mass_of_bits                                               ! mass of bergy bits (kg)
191            zdMbitsE = rn_bits_erosion_fraction * zdMe                               ! change in mass of bits (kg)
192            znMbits  = zMbits + zdMbitsE                                             ! add new bergy bits to mass (kg)
193            zLbits   = MIN( zL, zW, zT, 40._wp )                                     ! assume bergy bits are smallest dimension or 40 meters
194            zAbits   = ( zMbits / rn_rho_bergs ) / zLbits                            ! Effective bottom area (assuming T=Lbits)
195            zMbb     = MAX( 0.58_wp*(zdvo**0.8_wp)*(zSST+2._wp) /   &
196               &                              ( zLbits**0.2_wp ) , 0._wp ) * z1_rday ! Basal turbulent melting (for bits)
197            zMbb     = rn_rho_bergs * zAbits * zMbb                                  ! in kg/s
198            zdMbitsM = MIN( zMbb*zdt , znMbits )                                     ! bergy bits mass lost to melting (kg)
199            znMbits  = znMbits-zdMbitsM                                              ! remove mass lost to bergy bits melt
200            IF( zMnew == 0._wp ) THEN                                                ! if parent berg has completely melted then
201               zdMbitsM = zdMbitsM + znMbits                                         ! instantly melt all the bergy bits
202               znMbits  = 0._wp
203            ENDIF
204         ELSE                                                     ! No bergy bits
205            zAbits   = 0._wp
206            zdMbitsE = 0._wp
207            zdMbitsM = 0._wp
208            znMbits  = pt%mass_of_bits                             ! retain previous value incase non-zero
209         ENDIF
210
211         ! use tmask rather than tmask_i when dealing with icebergs
212         IF( tmask(ii,ij,1) /= 0._wp ) THEN    ! Add melting to the grid and field diagnostics
213            z1_e1e2    = r1_e1e2t(ii,ij) * this%mass_scaling
214            z1_dt_e1e2 = z1_dt * z1_e1e2
215            !
216            ! iceberg melt
217            !! the use of DDPDD function for the cumulative sum is needed for reproducibility
218            zmelt    = ( zdM - ( zdMbitsE - zdMbitsM ) ) * z1_dt   ! kg/s
219            CALL DDPDD( CMPLX( zmelt * z1_e1e2, 0.e0, dp ), cicb_melt(ii,ij) )
220            !
221            ! iceberg heat flux
222            !! the use of DDPDD function for the cumulative sum is needed for reproducibility
223            !! NB. The src_calving_hflx field is currently hardwired to zero in icb_stp, which means that the
224            !!     heat density of the icebergs is zero and the heat content flux to the ocean from iceberg
225            !!     melting is always zero. Leaving the term in the code until such a time as this is fixed. DS.
226            zheat_hcflux = zmelt * pt%heat_density       ! heat content flux : kg/s x J/kg = J/s
227            zheat_latent = - zmelt * rLfus               ! latent heat flux:  kg/s x J/kg = J/s
228            CALL DDPDD( CMPLX( ( zheat_hcflux + zheat_latent ) * z1_e1e2, 0.e0, dp ), cicb_hflx(ii,ij) )
229            !
230            ! diagnostics
231            CALL icb_dia_melt( ii, ij, zMnew, zheat_hcflux, zheat_latent, this%mass_scaling,       &
232               &                       zdM, zdMbitsE, zdMbitsM, zdMb, zdMe,   &
233               &                       zdMv, z1_dt_e1e2 )
234         ELSE
235            WRITE(numout,*) 'icb_thm: berg ',this%number(:),' appears to have grounded  at ',narea,ii,ij
236            CALL icb_utl_print_berg( this, kt )
237            WRITE(numout,*) 'msk=',tmask(ii,ij,1), e1e2t(ii,ij)
238            CALL ctl_stop('icb_thm', 'berg appears to have grounded!')
239         ENDIF
240
241         ! Rolling
242         zDn = rho_berg_1_oce * zTn       ! draught (keel depth)
243         IF( zDn > 0._wp .AND. MAX(zWn,zLn) < SQRT( 0.92*(zDn**2) + 58.32*zDn ) ) THEN
244            zT  = zTn
245            zTn = zWn
246            zWn = zT
247         ENDIF
248
249         ! Store the new state of iceberg (with L>W)
250         pt%mass         = zMnew
251         pt%mass_of_bits = znMbits
252         pt%thickness    = zTn
253         pt%width        = MIN( zWn , zLn )
254         pt%length       = MAX( zWn , zLn )
255
256         next=>this%next
257
258!!gm  add a test to avoid over melting ?
259!!pm  I agree, over melting could break conservation (more melt than calving)
260
261         IF( zMnew <= 0._wp ) THEN       ! Delete the berg if completely melted
262            CALL icb_utl_delete( first_berg, this )
263            !
264         ELSE                            ! Diagnose mass distribution on grid
265            z1_e1e2 = r1_e1e2t(ii,ij) * this%mass_scaling
266            CALL icb_dia_size( ii, ij, zWn, zLn, zAbits,   &
267               &               this%mass_scaling, zMnew, znMbits, z1_e1e2 )
268         ENDIF
269         !
270         this=>next
271         !
272      END DO
273      !
274      berg_grid%floating_melt = REAL(cicb_melt,dp)    ! kg/m2/s
275      berg_grid%calving_hflx  = REAL(cicb_hflx,dp)
276      !
277      ! now use melt and associated heat flux in ocean (or not)
278      !
279      IF(.NOT. ln_passive_mode ) THEN
280         emp (:,:) = emp (:,:) - berg_grid%floating_melt(:,:)
281         qns (:,:) = qns (:,:) + berg_grid%calving_hflx (:,:) 
282      ENDIF
283      !
284   END SUBROUTINE icb_thm
285
286   !!======================================================================
287END MODULE icbthm
Note: See TracBrowser for help on using the repository browser.