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

Last change on this file since 13226 was 13226, checked in by orioltp, 4 months ago

Merging dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation into the trunk

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