1 | MODULE 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 | !!---------------------------------------------------------------------- |
---|
37 | CONTAINS |
---|
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, 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 | berg_grid%calving_hflx(:,:) = 0._wp |
---|
71 | ! |
---|
72 | this => first_berg |
---|
73 | DO WHILE( ASSOCIATED(this) ) |
---|
74 | ! |
---|
75 | pt => this%current_point |
---|
76 | nknberg = this%number(1) |
---|
77 | CALL icb_utl_interp( pt%xi, pt%e1, pt%uo, pt%ui, pt%ua, pt%ssh_x, & |
---|
78 | & pt%yj, pt%e2, pt%vo, pt%vi, pt%va, pt%ssh_y, & |
---|
79 | & pt%sst, pt%cn, pt%hi, zff ) |
---|
80 | ! |
---|
81 | zSST = pt%sst |
---|
82 | zIC = MIN( 1._wp, pt%cn + rn_sicn_shift ) ! Shift sea-ice concentration !!gm ??? |
---|
83 | zM = pt%mass |
---|
84 | zT = pt%thickness ! total thickness |
---|
85 | ! D = (rn_rho_bergs/pp_rho_seawater)*zT ! draught (keel depth) |
---|
86 | ! F = zT - D ! freeboard |
---|
87 | zW = pt%width |
---|
88 | zL = pt%length |
---|
89 | zxi = pt%xi ! position in (i,j) referential |
---|
90 | zyj = pt%yj |
---|
91 | ii = INT( zxi + 0.5 ) ! T-cell of the berg |
---|
92 | ii = mi1( ii ) |
---|
93 | ij = INT( zyj + 0.5 ) |
---|
94 | ij = mj1( ij ) |
---|
95 | zVol = zT * zW * zL |
---|
96 | |
---|
97 | ! Environment |
---|
98 | zdvo = SQRT( (pt%uvel-pt%uo)**2 + (pt%vvel-pt%vo)**2 ) |
---|
99 | zdva = SQRT( (pt%ua -pt%uo)**2 + (pt%va -pt%vo)**2 ) |
---|
100 | zSs = 1.5_wp * SQRT( zdva ) + 0.1_wp * zdva ! Sea state (eqn M.A9) |
---|
101 | |
---|
102 | ! Melt rates in m/s (i.e. division by rday) |
---|
103 | zMv = MAX( 7.62d-3*zSST+1.29d-3*(zSST**2) , 0._wp ) * z1_rday ! Buoyant convection at sides (eqn M.A10) |
---|
104 | 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 ) |
---|
105 | zMe = MAX( z1_12*(zSST+2.)*zSs*(1._wp+COS(rpi*(zIC**3))) , 0._wp ) * z1_rday ! Wave erosion (eqn M.A8 ) |
---|
106 | |
---|
107 | IF( ln_operator_splitting ) THEN ! Operator split update of volume/mass |
---|
108 | zTn = MAX( zT - zMb*zdt , 0._wp ) ! new total thickness (m) |
---|
109 | znVol = zTn * zW * zL ! new volume (m^3) |
---|
110 | zMnew1 = ( znVol / zVol ) * zM ! new mass (kg) |
---|
111 | zdMb = zM - zMnew1 ! mass lost to basal melting (>0) (kg) |
---|
112 | ! |
---|
113 | zLn = MAX( zL - zMv*zdt , 0._wp ) ! new length (m) |
---|
114 | zWn = MAX( zW - zMv*zdt , 0._wp ) ! new width (m) |
---|
115 | znVol = zTn * zWn * zLn ! new volume (m^3) |
---|
116 | zMnew2 = ( znVol / zVol ) * zM ! new mass (kg) |
---|
117 | zdMv = zMnew1 - zMnew2 ! mass lost to buoyant convection (>0) (kg) |
---|
118 | ! |
---|
119 | zLn = MAX( zLn - zMe*zdt , 0._wp ) ! new length (m) |
---|
120 | zWn = MAX( zWn - zMe*zdt , 0._wp ) ! new width (m) |
---|
121 | znVol = zTn * zWn * zLn ! new volume (m^3) |
---|
122 | zMnew = ( znVol / zVol ) * zM ! new mass (kg) |
---|
123 | zdMe = zMnew2 - zMnew ! mass lost to erosion (>0) (kg) |
---|
124 | zdM = zM - zMnew ! mass lost to all erosion and melting (>0) (kg) |
---|
125 | ! |
---|
126 | ELSE ! Update dimensions of berg |
---|
127 | zLn = MAX( zL -(zMv+zMe)*zdt ,0._wp ) ! (m) |
---|
128 | zWn = MAX( zW -(zMv+zMe)*zdt ,0._wp ) ! (m) |
---|
129 | zTn = MAX( zT - zMb *zdt ,0._wp ) ! (m) |
---|
130 | ! Update volume and mass of berg |
---|
131 | znVol = zTn*zWn*zLn ! (m^3) |
---|
132 | zMnew = (znVol/zVol)*zM ! (kg) |
---|
133 | zdM = zM - zMnew ! (kg) |
---|
134 | zdMb = (zM/zVol) * (zW* zL ) *zMb*zdt ! approx. mass loss to basal melting (kg) |
---|
135 | zdMe = (zM/zVol) * (zT*(zW+zL)) *zMe*zdt ! approx. mass lost to erosion (kg) |
---|
136 | zdMv = (zM/zVol) * (zT*(zW+zL)) *zMv*zdt ! approx. mass loss to buoyant convection (kg) |
---|
137 | ENDIF |
---|
138 | |
---|
139 | IF( rn_bits_erosion_fraction > 0._wp ) THEN ! Bergy bits |
---|
140 | ! |
---|
141 | zMbits = pt%mass_of_bits ! mass of bergy bits (kg) |
---|
142 | zdMbitsE = rn_bits_erosion_fraction * zdMe ! change in mass of bits (kg) |
---|
143 | znMbits = zMbits + zdMbitsE ! add new bergy bits to mass (kg) |
---|
144 | zLbits = MIN( zL, zW, zT, 40._wp ) ! assume bergy bits are smallest dimension or 40 meters |
---|
145 | zAbits = ( zMbits / rn_rho_bergs ) / zLbits ! Effective bottom area (assuming T=Lbits) |
---|
146 | zMbb = MAX( 0.58_wp*(zdvo**0.8_wp)*(zSST+2._wp) / & |
---|
147 | & ( zLbits**0.2_wp ) , 0._wp ) * z1_rday ! Basal turbulent melting (for bits) |
---|
148 | zMbb = rn_rho_bergs * zAbits * zMbb ! in kg/s |
---|
149 | zdMbitsM = MIN( zMbb*zdt , znMbits ) ! bergy bits mass lost to melting (kg) |
---|
150 | znMbits = znMbits-zdMbitsM ! remove mass lost to bergy bits melt |
---|
151 | IF( zMnew == 0._wp ) THEN ! if parent berg has completely melted then |
---|
152 | zdMbitsM = zdMbitsM + znMbits ! instantly melt all the bergy bits |
---|
153 | znMbits = 0._wp |
---|
154 | ENDIF |
---|
155 | ELSE ! No bergy bits |
---|
156 | zAbits = 0._wp |
---|
157 | zdMbitsE = 0._wp |
---|
158 | zdMbitsM = 0._wp |
---|
159 | znMbits = pt%mass_of_bits ! retain previous value incase non-zero |
---|
160 | ENDIF |
---|
161 | |
---|
162 | ! use tmask rather than tmask_i when dealing with icebergs |
---|
163 | IF( tmask(ii,ij,1) /= 0._wp ) THEN ! Add melting to the grid and field diagnostics |
---|
164 | z1_e1e2 = r1_e1e2t(ii,ij) * this%mass_scaling |
---|
165 | z1_dt_e1e2 = z1_dt * z1_e1e2 |
---|
166 | zmelt = ( zdM - ( zdMbitsE - zdMbitsM ) ) * z1_dt ! kg/s |
---|
167 | berg_grid%floating_melt(ii,ij) = berg_grid%floating_melt(ii,ij) + zmelt * z1_e1e2 ! kg/m2/s |
---|
168 | zheat = zmelt * pt%heat_density ! kg/s x J/kg = J/s |
---|
169 | berg_grid%calving_hflx (ii,ij) = berg_grid%calving_hflx (ii,ij) + zheat * z1_e1e2 ! W/m2 |
---|
170 | CALL icb_dia_melt( ii, ij, zMnew, zheat, this%mass_scaling, & |
---|
171 | & zdM, zdMbitsE, zdMbitsM, zdMb, zdMe, & |
---|
172 | & zdMv, z1_dt_e1e2 ) |
---|
173 | ELSE |
---|
174 | WRITE(numout,*) 'icb_thm: berg ',this%number(:),' appears to have grounded at ',narea,ii,ij |
---|
175 | CALL icb_utl_print_berg( this, kt ) |
---|
176 | WRITE(numout,*) 'msk=',tmask(ii,ij,1), e1e2t(ii,ij) |
---|
177 | CALL ctl_stop('icb_thm', 'berg appears to have grounded!') |
---|
178 | ENDIF |
---|
179 | |
---|
180 | ! Rolling |
---|
181 | zDn = ( rn_rho_bergs / pp_rho_seawater ) * zTn ! draught (keel depth) |
---|
182 | IF( zDn > 0._wp .AND. MAX(zWn,zLn) < SQRT( 0.92*(zDn**2) + 58.32*zDn ) ) THEN |
---|
183 | zT = zTn |
---|
184 | zTn = zWn |
---|
185 | zWn = zT |
---|
186 | ENDIF |
---|
187 | |
---|
188 | ! Store the new state of iceberg (with L>W) |
---|
189 | pt%mass = zMnew |
---|
190 | pt%mass_of_bits = znMbits |
---|
191 | pt%thickness = zTn |
---|
192 | pt%width = MIN( zWn , zLn ) |
---|
193 | pt%length = MAX( zWn , zLn ) |
---|
194 | |
---|
195 | next=>this%next |
---|
196 | |
---|
197 | !!gm add a test to avoid over melting ? |
---|
198 | |
---|
199 | IF( zMnew <= 0._wp ) THEN ! Delete the berg if completely melted |
---|
200 | CALL icb_utl_delete( first_berg, this ) |
---|
201 | ! |
---|
202 | ELSE ! Diagnose mass distribution on grid |
---|
203 | z1_e1e2 = r1_e1e2t(ii,ij) * this%mass_scaling |
---|
204 | CALL icb_dia_size( ii, ij, zWn, zLn, zAbits, & |
---|
205 | & this%mass_scaling, zMnew, znMbits, z1_e1e2 ) |
---|
206 | ENDIF |
---|
207 | ! |
---|
208 | this=>next |
---|
209 | ! |
---|
210 | END DO |
---|
211 | |
---|
212 | ! now use melt and associated heat flux in ocean (or not) |
---|
213 | ! |
---|
214 | IF(.NOT. ln_passive_mode ) THEN |
---|
215 | emp (:,:) = emp (:,:) - berg_grid%floating_melt(:,:) |
---|
216 | !! qns (:,:) = qns (:,:) + berg_grid%calving_hflx (:,:) !!gm heat flux not yet properly coded ==>> need it, SOLVE that! |
---|
217 | ENDIF |
---|
218 | ! |
---|
219 | END SUBROUTINE icb_thm |
---|
220 | |
---|
221 | !!====================================================================== |
---|
222 | END MODULE icbthm |
---|