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

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

Last change on this file since 10068 was 10068, checked in by nicolasmartin, 6 years ago

First part of modifications to have a common default header : fix typos and SVN keywords properties

  • Property svn:keywords set to Id
File size: 12.1 KB
RevLine 
[3614]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
[9190]32   !!----------------------------------------------------------------------
[9598]33   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[5215]34   !! $Id$
[10068]35   !! Software governed by the CeCILL license (see ./LICENSE)
[9190]36   !!----------------------------------------------------------------------
[3614]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
[9990]52      REAL(wp) ::   zMnew, zMnew1, zMnew2, zheat_hcflux, zheat_latent, z1_12
[3614]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
[9940]60      z1_12   = 1._wp / 12._wp
61      zdt     = berg_dt
62      z1_dt   = 1._wp / zdt
[9190]63      !
[3614]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
[9990]70      ! calving_hflx re-used here as temporary workspace for the heat flux associated with melting
[3614]71      berg_grid%calving_hflx(:,:)  = 0._wp
[9190]72      !
[3614]73      this => first_berg
[9190]74      DO WHILE( ASSOCIATED(this) )
[3614]75         !
76         pt => this%current_point
77         nknberg = this%number(1)
[9190]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 )
[3614]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 )
[9940]101         zSs  = 1.5_wp * SQRT( zdva ) + 0.1_wp * zdva                ! Sea state      (eqn M.A9)
[3614]102
103         ! Melt rates in m/s (i.e. division by rday)
[9940]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 )
[3614]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)
[9190]110            znVol  = zTn * zW * zL                       ! new volume (m^3)
[9940]111            zMnew1 = ( znVol / zVol ) * zM               ! new mass (kg)
[3614]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)
[9190]116            znVol  = zTn * zWn * zLn                     ! new volume (m^3)
[9940]117            zMnew2 = ( znVol / zVol ) * zM               ! new mass (kg)
[3614]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)
[9190]122            znVol  = zTn * zWn * zLn                     ! new volume (m^3)
123            zMnew  = ( znVol / zVol ) * zM               ! new mass (kg)
[3614]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
[9190]128            zLn = MAX( zL -(zMv+zMe)*zdt ,0._wp )        ! (m)
129            zWn = MAX( zW -(zMv+zMe)*zdt ,0._wp )        ! (m)
[3614]130            zTn = MAX( zT - zMb    *zdt ,0._wp )         ! (m)
131            ! Update volume and mass of berg
[9190]132            znVol = zTn*zWn*zLn                          ! (m^3)
133            zMnew = (znVol/zVol)*zM                      ! (kg)
[3614]134            zdM   = zM - zMnew                           ! (kg)
[9190]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)
[3614]138         ENDIF
139
[9190]140         IF( rn_bits_erosion_fraction > 0._wp ) THEN     ! Bergy bits
[3614]141            !
142            zMbits   = pt%mass_of_bits                                               ! mass of bergy bits (kg)
[9190]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)
[9940]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)
[9190]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
[3614]152            IF( zMnew == 0._wp ) THEN                                                ! if parent berg has completely melted then
[9190]153               zdMbitsM = zdMbitsM + znMbits                                         ! instantly melt all the bergy bits
[3614]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
[5836]165            z1_e1e2    = r1_e1e2t(ii,ij) * this%mass_scaling
[3614]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
[9990]169            !! NB. The src_calving_hflx field is currently hardwired to zero in icb_stp, which means that the
170            !!     heat density of the icebergs is zero and the heat content flux to the ocean from iceberg
171            !!     melting is always zero. Leaving the term in the code until such a time as this is fixed. DS.
172            zheat_hcflux = zmelt * pt%heat_density       ! heat content flux : kg/s x J/kg = J/s
173            zheat_latent = - zmelt * rLfus               ! latent heat flux:  kg/s x J/kg = J/s
174            berg_grid%calving_hflx (ii,ij) = berg_grid%calving_hflx (ii,ij) + ( zheat_hcflux + zheat_latent ) * z1_e1e2    ! W/m2
175            CALL icb_dia_melt( ii, ij, zMnew, zheat_hcflux, zheat_latent, this%mass_scaling,       &
[9190]176               &                       zdM, zdMbitsE, zdMbitsM, zdMb, zdMe,   &
177               &                       zdMv, z1_dt_e1e2 )
[3614]178         ELSE
179            WRITE(numout,*) 'icb_thm: berg ',this%number(:),' appears to have grounded  at ',narea,ii,ij
180            CALL icb_utl_print_berg( this, kt )
181            WRITE(numout,*) 'msk=',tmask(ii,ij,1), e1e2t(ii,ij)
182            CALL ctl_stop('icb_thm', 'berg appears to have grounded!')
183         ENDIF
184
185         ! Rolling
186         zDn = ( rn_rho_bergs / pp_rho_seawater ) * zTn       ! draught (keel depth)
187         IF( zDn > 0._wp .AND. MAX(zWn,zLn) < SQRT( 0.92*(zDn**2) + 58.32*zDn ) ) THEN
188            zT  = zTn
189            zTn = zWn
190            zWn = zT
[9190]191         ENDIF
[3614]192
193         ! Store the new state of iceberg (with L>W)
194         pt%mass         = zMnew
195         pt%mass_of_bits = znMbits
196         pt%thickness    = zTn
[9190]197         pt%width        = MIN( zWn , zLn )
198         pt%length       = MAX( zWn , zLn )
[3614]199
200         next=>this%next
201
202!!gm  add a test to avoid over melting ?
203
204         IF( zMnew <= 0._wp ) THEN       ! Delete the berg if completely melted
205            CALL icb_utl_delete( first_berg, this )
206            !
207         ELSE                            ! Diagnose mass distribution on grid
[5836]208            z1_e1e2 = r1_e1e2t(ii,ij) * this%mass_scaling
[3614]209            CALL icb_dia_size( ii, ij, zWn, zLn, zAbits,   &
[9190]210               &               this%mass_scaling, zMnew, znMbits, z1_e1e2 )
[3614]211         ENDIF
212         !
213         this=>next
214         !
215      END DO
[9190]216
[3614]217      ! now use melt and associated heat flux in ocean (or not)
218      !
219      IF(.NOT. ln_passive_mode ) THEN
220         emp (:,:) = emp (:,:) - berg_grid%floating_melt(:,:)
[9990]221         qns (:,:) = qns (:,:) + berg_grid%calving_hflx (:,:) 
[3614]222      ENDIF
223      !
224   END SUBROUTINE icb_thm
225
226   !!======================================================================
227END MODULE icbthm
Note: See TracBrowser for help on using the repository browser.