source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/ICB/icbthm.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 13 months ago

The Dr Hook changes from my perl code.

File size: 11.6 KB
Line 
1MODULE icbthm
2
3   !!======================================================================
4   !!                       ***  MODULE  icbthm  ***
5   !! Icebergs:  thermodynamics routines for icebergs
6   !!======================================================================
7   !! History : 3.3.1 !  2010-01  (Martin&Adcroft) Original code
8   !!            -    !  2011-03  (Madec)          Part conversion to NEMO form
9   !!            -    !                            Removal of mapping from another grid
10   !!            -    !  2011-04  (Alderson)       Split into separate modules
11   !!            -    !  2011-05  (Alderson)       Use tmask instead of tmask_i
12   !!----------------------------------------------------------------------
13   !!----------------------------------------------------------------------
14   !!   icb_thm : initialise
15   !!             reference for equations - M = Martin + Adcroft, OM 34, 2010
16   !!----------------------------------------------------------------------
17   USE par_oce        ! NEMO parameters
18   USE dom_oce        ! NEMO domain
19   USE in_out_manager ! NEMO IO routines, numout in particular
20   USE iom
21   USE lib_mpp        ! NEMO MPI routines, ctl_stop in particular
22   USE phycst         ! NEMO physical constants
23   USE sbc_oce
24
25   USE icb_oce        ! define iceberg arrays
26   USE icbutl         ! iceberg utility routines
27   USE icbdia         ! iceberg budget routines
28
29   USE yomhook, ONLY: lhook, dr_hook
30   USE parkind1, ONLY: jprb, jpim
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   icb_thm ! routine called in icbstp.F90 module
36
37   !! $Id$
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
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      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
59      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
60      REAL(KIND=jprb)               :: zhook_handle
61
62      CHARACTER(LEN=*), PARAMETER :: RoutineName='ICB_THM'
63
64      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
65
66      !!----------------------------------------------------------------------
67      !
68      z1_rday = 1._wp / rday
69     
70      ! we're either going to ignore berg fresh water melt flux and associated heat
71      ! or we pass it into the ocean, so at this point we set them both to zero,
72      ! accumulate the contributions to them from each iceberg in the while loop following
73      ! and then pass them (or not) to the ocean
74      !
75      berg_grid%floating_melt(:,:) = 0._wp
76      berg_grid%calving_hflx(:,:)  = 0._wp
77   
78      this => first_berg
79      DO WHILE( associated(this) )
80         !
81         pt => this%current_point
82         nknberg = this%number(1)
83         CALL icb_utl_interp( pt%xi, pt%e1, pt%uo, pt%ui, pt%ua, pt%ssh_x, &
84         &                    pt%yj, pt%e2, pt%vo, pt%vi, pt%va, pt%ssh_y, &
85         &                 pt%sst, pt%cn, pt%hi, zff )
86         !
87         zSST = pt%sst
88         zIC  = MIN( 1._wp, pt%cn + rn_sicn_shift )     ! Shift sea-ice concentration       !!gm ???
89         zM   = pt%mass
90         zT   = pt%thickness                               ! total thickness
91       ! D   = (rn_rho_bergs/pp_rho_seawater)*zT ! draught (keel depth)
92       ! F   = zT - D ! freeboard
93         zW   = pt%width
94         zL   = pt%length
95         zxi  = pt%xi                                      ! position in (i,j) referential
96         zyj  = pt%yj
97         ii  = INT( zxi + 0.5 )                            ! T-cell of the berg
98         ii  = mi1( ii )
99         ij  = INT( zyj + 0.5 )             
100         ij  = mj1( ij )
101         zVol = zT * zW * zL
102         zdt = berg_dt   ;   z1_dt = 1._wp / zdt
103
104         ! Environment
105         zdvo = SQRT( (pt%uvel-pt%uo)**2 + (pt%vvel-pt%vo)**2 )
106         zdva = SQRT( (pt%ua  -pt%uo)**2 + (pt%va  -pt%vo)**2 )
107         zSs  = 1.5 * SQRT( zdva ) + 0.1 * zdva                ! Sea state      (eqn M.A9)
108
109         ! Melt rates in m/s (i.e. division by rday)
110         zMv = MAX( 7.62e-3*zSST+1.29e-3*(zSST**2)            , 0._wp ) * z1_rday   ! Buoyant convection at sides (eqn M.A10)
111         zMb = MAX( 0.58*(zdvo**0.8)*(zSST+4.0)/(zL**0.2)      , 0._wp ) * z1_rday   ! Basal turbulent melting     (eqn M.A7 )
112         zMe = MAX( 1./12.*(zSST+2.)*zSs*(1+cos(rpi*(zIC**3))) , 0._wp ) * z1_rday   ! Wave erosion                (eqn M.A8 )
113
114         IF( ln_operator_splitting ) THEN      ! Operator split update of volume/mass
115            zTn    = MAX( zT - zMb*zdt , 0._wp )         ! new total thickness (m)
116            znVol  = zTn * zW * zL                        ! new volume (m^3)
117            zMnew1 = (znVol/zVol) * zM                    ! new mass (kg)
118            zdMb   = zM - zMnew1                         ! mass lost to basal melting (>0) (kg)
119            !
120            zLn    = MAX( zL - zMv*zdt , 0._wp )         ! new length (m)
121            zWn    = MAX( zW - zMv*zdt , 0._wp )         ! new width (m)
122            znVol  = zTn * zWn * zLn                      ! new volume (m^3)
123            zMnew2 = (znVol/zVol) * zM                    ! new mass (kg)
124            zdMv   = zMnew1 - zMnew2                     ! mass lost to buoyant convection (>0) (kg)
125            !
126            zLn    = MAX( zLn - zMe*zdt , 0._wp )        ! new length (m)
127            zWn    = MAX( zWn - zMe*zdt , 0._wp )        ! new width (m)
128            znVol  = zTn * zWn * zLn                      ! new volume (m^3)
129            zMnew  = ( znVol / zVol ) * zM                ! new mass (kg)
130            zdMe   = zMnew2 - zMnew                      ! mass lost to erosion (>0) (kg)
131            zdM    = zM - zMnew                          ! mass lost to all erosion and melting (>0) (kg)
132            !
133         ELSE                                         ! Update dimensions of berg
134            zLn = MAX( zL -(zMv+zMe)*zdt ,0._wp )         ! (m)
135            zWn = MAX( zW -(zMv+zMe)*zdt ,0._wp )         ! (m)
136            zTn = MAX( zT - zMb    *zdt ,0._wp )         ! (m)
137            ! Update volume and mass of berg
138            znVol = zTn*zWn*zLn                           ! (m^3)
139            zMnew = (znVol/zVol)*zM                       ! (kg)
140            zdM   = zM - zMnew                           ! (kg)
141            zdMb = (zM/zVol) * (zW*   zL ) *zMb*zdt         ! approx. mass loss to basal melting (kg)
142            zdMe = (zM/zVol) * (zT*(zW+zL)) *zMe*zdt         ! approx. mass lost to erosion (kg)
143            zdMv = (zM/zVol) * (zT*(zW+zL)) *zMv*zdt         ! approx. mass loss to buoyant convection (kg)
144         ENDIF
145
146         IF( rn_bits_erosion_fraction > 0._wp ) THEN      ! Bergy bits
147            !
148            zMbits   = pt%mass_of_bits                                               ! mass of bergy bits (kg)
149            zdMbitsE = rn_bits_erosion_fraction * zdMe                        ! change in mass of bits (kg)
150            znMbits  = zMbits + zdMbitsE                                               ! add new bergy bits to mass (kg)
151            zLbits   = MIN( zL, zW, zT, 40._wp )                                        ! assume bergy bits are smallest dimension or 40 meters
152            zAbits   = ( zMbits / rn_rho_bergs ) / zLbits                           ! Effective bottom area (assuming T=Lbits)
153            zMbb     = MAX( 0.58*(zdvo**0.8)*(zSST+2.0)/(zLbits**0.2), 0.) * z1_rday    ! Basal turbulent melting (for bits)
154            zMbb     = rn_rho_bergs * zAbits * zMbb                                 ! in kg/s
155            zdMbitsM = MIN( zMbb*zdt , znMbits )                                       ! bergy bits mass lost to melting (kg)
156            znMbits  = znMbits-zdMbitsM                                                ! remove mass lost to bergy bits melt
157            IF( zMnew == 0._wp ) THEN                                                ! if parent berg has completely melted then
158               zdMbitsM = zdMbitsM + znMbits                                           ! instantly melt all the bergy bits
159               znMbits  = 0._wp
160            ENDIF
161         ELSE                                                     ! No bergy bits
162            zAbits   = 0._wp
163            zdMbitsE = 0._wp
164            zdMbitsM = 0._wp
165            znMbits  = pt%mass_of_bits                             ! retain previous value incase non-zero
166         ENDIF
167
168         ! use tmask rather than tmask_i when dealing with icebergs
169         IF( tmask(ii,ij,1) /= 0._wp ) THEN    ! Add melting to the grid and field diagnostics
170            z1_e1e2    = 1._wp / e1e2t(ii,ij) * this%mass_scaling
171            z1_dt_e1e2 = z1_dt * z1_e1e2
172            zmelt    = ( zdM - ( zdMbitsE - zdMbitsM ) ) * z1_dt   ! kg/s
173            berg_grid%floating_melt(ii,ij) = berg_grid%floating_melt(ii,ij) + zmelt    * z1_e1e2    ! kg/m2/s
174!            zheat = zmelt * pt%heat_density              ! kg/s x J/kg = J/s
175            zheat = zmelt * lfus                           !rma kg/s x J/kg (latent heat of fusion) = J/s
176            berg_grid%calving_hflx (ii,ij) = berg_grid%calving_hflx (ii,ij) + zheat    * z1_e1e2    ! W/m2
177            CALL icb_dia_melt( ii, ij, zMnew, zheat, this%mass_scaling,       &
178            &                          zdM, zdMbitsE, zdMbitsM, zdMb, zdMe,   &
179            &                          zdMv, z1_dt_e1e2 )
180         ELSE
181            WRITE(numout,*) 'icb_thm: berg ',this%number(:),' appears to have grounded  at ',narea,ii,ij
182            CALL icb_utl_print_berg( this, kt )
183            WRITE(numout,*) 'msk=',tmask(ii,ij,1), e1e2t(ii,ij)
184            CALL ctl_stop('icb_thm', 'berg appears to have grounded!')
185         ENDIF
186
187         ! Rolling
188         zDn = ( rn_rho_bergs / pp_rho_seawater ) * zTn       ! draught (keel depth)
189         IF( zDn > 0._wp .AND. MAX(zWn,zLn) < SQRT( 0.92*(zDn**2) + 58.32*zDn ) ) THEN
190            zT  = zTn
191            zTn = zWn
192            zWn = zT
193         endif
194
195         ! Store the new state of iceberg (with L>W)
196         pt%mass         = zMnew
197         pt%mass_of_bits = znMbits
198         pt%thickness    = zTn
199         pt%width        = min(zWn,zLn)
200         pt%length       = max(zWn,zLn)
201
202         next=>this%next
203
204!!gm  add a test to avoid over melting ?
205
206         IF( zMnew <= 0._wp ) THEN       ! Delete the berg if completely melted
207            CALL icb_utl_delete( first_berg, this )
208            !
209         ELSE                            ! Diagnose mass distribution on grid
210            z1_e1e2 = 1._wp / e1e2t(ii,ij) * this%mass_scaling
211            CALL icb_dia_size( ii, ij, zWn, zLn, zAbits,   &
212            &                  this%mass_scaling, zMnew, znMbits, z1_e1e2)
213         ENDIF
214         !
215         this=>next
216         !
217      END DO
218     
219      ! now use melt and associated heat flux in ocean (or not)
220      !
221      IF(.NOT. ln_passive_mode ) THEN
222         emp (:,:) = emp (:,:) - berg_grid%floating_melt(:,:)
223         qns (:,:) = qns (:,:) - berg_grid%calving_hflx (:,:) 
224      ENDIF
225      !
226      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
227   END SUBROUTINE icb_thm
228
229   !!======================================================================
230END MODULE icbthm
Note: See TracBrowser for help on using the repository browser.