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/icebergs_ocean_heat_fluxes/src/OCE/ICB – NEMO

source: NEMO/branches/UKMO/icebergs_ocean_heat_fluxes/src/OCE/ICB/icbthm.F90 @ 9959

Last change on this file since 9959 was 9959, checked in by davestorkey, 6 years ago

UKMO icebergs_ocean_heat_fluxes branch : science changes

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