- Timestamp:
- 2018-01-06T15:18:23+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icbthm.F90
r5836 r9190 1 1 MODULE icbthm 2 3 2 !!====================================================================== 4 3 !! *** MODULE icbthm *** … … 31 30 PUBLIC icb_thm ! routine called in icbstp.F90 module 32 31 32 !!---------------------------------------------------------------------- 33 !! NEMO/OPA 3.3 , NEMO Consortium (2011) 33 34 !! $Id$ 35 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 36 !!---------------------------------------------------------------------- 34 37 CONTAINS 35 38 … … 55 58 ! 56 59 z1_rday = 1._wp / rday 57 60 ! 58 61 ! we're either going to ignore berg fresh water melt flux and associated heat 59 62 ! or we pass it into the ocean, so at this point we set them both to zero, … … 63 66 berg_grid%floating_melt(:,:) = 0._wp 64 67 berg_grid%calving_hflx(:,:) = 0._wp 65 68 ! 66 69 this => first_berg 67 DO WHILE( associated(this) )70 DO WHILE( ASSOCIATED(this) ) 68 71 ! 69 72 pt => this%current_point 70 73 nknberg = this%number(1) 71 CALL icb_utl_interp( pt%xi, pt%e1, pt%uo, pt%ui, pt%ua, pt%ssh_x, &72 & pt%yj, pt%e2, pt%vo, pt%vi, pt%va, pt%ssh_y,&73 & pt%sst, pt%cn, pt%hi, zff )74 CALL icb_utl_interp( pt%xi, pt%e1, pt%uo, pt%ui, pt%ua, pt%ssh_x, & 75 & pt%yj, pt%e2, pt%vo, pt%vi, pt%va, pt%ssh_y, & 76 & pt%sst, pt%cn, pt%hi, zff ) 74 77 ! 75 78 zSST = pt%sst … … 98 101 zMv = MAX( 7.62e-3*zSST+1.29e-3*(zSST**2) , 0._wp ) * z1_rday ! Buoyant convection at sides (eqn M.A10) 99 102 zMb = MAX( 0.58*(zdvo**0.8)*(zSST+4.0)/(zL**0.2) , 0._wp ) * z1_rday ! Basal turbulent melting (eqn M.A7 ) 100 zMe = MAX( 1./12.*(zSST+2.)*zSs*(1+ cos(rpi*(zIC**3))) , 0._wp ) * z1_rday ! Wave erosion (eqn M.A8 )103 zMe = MAX( 1./12.*(zSST+2.)*zSs*(1+COS(rpi*(zIC**3))) , 0._wp ) * z1_rday ! Wave erosion (eqn M.A8 ) 101 104 102 105 IF( ln_operator_splitting ) THEN ! Operator split update of volume/mass 103 106 zTn = MAX( zT - zMb*zdt , 0._wp ) ! new total thickness (m) 104 znVol = zTn * zW * zL 105 zMnew1 = (znVol/zVol) * zM 107 znVol = zTn * zW * zL ! new volume (m^3) 108 zMnew1 = (znVol/zVol) * zM ! new mass (kg) 106 109 zdMb = zM - zMnew1 ! mass lost to basal melting (>0) (kg) 107 110 ! 108 111 zLn = MAX( zL - zMv*zdt , 0._wp ) ! new length (m) 109 112 zWn = MAX( zW - zMv*zdt , 0._wp ) ! new width (m) 110 znVol = zTn * zWn * zLn 111 zMnew2 = (znVol/zVol) * zM 113 znVol = zTn * zWn * zLn ! new volume (m^3) 114 zMnew2 = (znVol/zVol) * zM ! new mass (kg) 112 115 zdMv = zMnew1 - zMnew2 ! mass lost to buoyant convection (>0) (kg) 113 116 ! 114 117 zLn = MAX( zLn - zMe*zdt , 0._wp ) ! new length (m) 115 118 zWn = MAX( zWn - zMe*zdt , 0._wp ) ! new width (m) 116 znVol = zTn * zWn * zLn 117 zMnew = ( znVol / zVol ) * zM 119 znVol = zTn * zWn * zLn ! new volume (m^3) 120 zMnew = ( znVol / zVol ) * zM ! new mass (kg) 118 121 zdMe = zMnew2 - zMnew ! mass lost to erosion (>0) (kg) 119 122 zdM = zM - zMnew ! mass lost to all erosion and melting (>0) (kg) 120 123 ! 121 124 ELSE ! Update dimensions of berg 122 zLn = MAX( zL -(zMv+zMe)*zdt ,0._wp ) 123 zWn = MAX( zW -(zMv+zMe)*zdt ,0._wp ) 125 zLn = MAX( zL -(zMv+zMe)*zdt ,0._wp ) ! (m) 126 zWn = MAX( zW -(zMv+zMe)*zdt ,0._wp ) ! (m) 124 127 zTn = MAX( zT - zMb *zdt ,0._wp ) ! (m) 125 128 ! Update volume and mass of berg 126 znVol = zTn*zWn*zLn 127 zMnew = (znVol/zVol)*zM 129 znVol = zTn*zWn*zLn ! (m^3) 130 zMnew = (znVol/zVol)*zM ! (kg) 128 131 zdM = zM - zMnew ! (kg) 129 zdMb = (zM/zVol) * (zW* zL ) *zMb*zdt 130 zdMe = (zM/zVol) * (zT*(zW+zL)) *zMe*zdt 131 zdMv = (zM/zVol) * (zT*(zW+zL)) *zMv*zdt 132 ENDIF 133 134 IF( rn_bits_erosion_fraction > 0._wp ) THEN 132 zdMb = (zM/zVol) * (zW* zL ) *zMb*zdt ! approx. mass loss to basal melting (kg) 133 zdMe = (zM/zVol) * (zT*(zW+zL)) *zMe*zdt ! approx. mass lost to erosion (kg) 134 zdMv = (zM/zVol) * (zT*(zW+zL)) *zMv*zdt ! approx. mass loss to buoyant convection (kg) 135 ENDIF 136 137 IF( rn_bits_erosion_fraction > 0._wp ) THEN ! Bergy bits 135 138 ! 136 139 zMbits = pt%mass_of_bits ! mass of bergy bits (kg) 137 zdMbitsE = rn_bits_erosion_fraction * zdMe ! change in mass of bits (kg)138 znMbits = zMbits + zdMbitsE 139 zLbits = MIN( zL, zW, zT, 40._wp ) 140 zAbits = ( zMbits / rn_rho_bergs ) / zLbits ! Effective bottom area (assuming T=Lbits)141 zMbb = MAX( 0.58*(zdvo**0.8)*(zSST+2.0)/(zLbits**0.2), 0.) * z1_rday 142 zMbb = rn_rho_bergs * zAbits * zMbb ! in kg/s143 zdMbitsM = MIN( zMbb*zdt , znMbits ) 144 znMbits = znMbits-zdMbitsM 140 zdMbitsE = rn_bits_erosion_fraction * zdMe ! change in mass of bits (kg) 141 znMbits = zMbits + zdMbitsE ! add new bergy bits to mass (kg) 142 zLbits = MIN( zL, zW, zT, 40._wp ) ! assume bergy bits are smallest dimension or 40 meters 143 zAbits = ( zMbits / rn_rho_bergs ) / zLbits ! Effective bottom area (assuming T=Lbits) 144 zMbb = MAX( 0.58*(zdvo**0.8)*(zSST+2.0)/(zLbits**0.2), 0.) * z1_rday ! Basal turbulent melting (for bits) 145 zMbb = rn_rho_bergs * zAbits * zMbb ! in kg/s 146 zdMbitsM = MIN( zMbb*zdt , znMbits ) ! bergy bits mass lost to melting (kg) 147 znMbits = znMbits-zdMbitsM ! remove mass lost to bergy bits melt 145 148 IF( zMnew == 0._wp ) THEN ! if parent berg has completely melted then 146 zdMbitsM = zdMbitsM + znMbits 149 zdMbitsM = zdMbitsM + znMbits ! instantly melt all the bergy bits 147 150 znMbits = 0._wp 148 151 ENDIF … … 163 166 berg_grid%calving_hflx (ii,ij) = berg_grid%calving_hflx (ii,ij) + zheat * z1_e1e2 ! W/m2 164 167 CALL icb_dia_melt( ii, ij, zMnew, zheat, this%mass_scaling, & 165 &zdM, zdMbitsE, zdMbitsM, zdMb, zdMe, &166 &zdMv, z1_dt_e1e2 )168 & zdM, zdMbitsE, zdMbitsM, zdMb, zdMe, & 169 & zdMv, z1_dt_e1e2 ) 167 170 ELSE 168 171 WRITE(numout,*) 'icb_thm: berg ',this%number(:),' appears to have grounded at ',narea,ii,ij … … 178 181 zTn = zWn 179 182 zWn = zT 180 endif183 ENDIF 181 184 182 185 ! Store the new state of iceberg (with L>W) … … 184 187 pt%mass_of_bits = znMbits 185 188 pt%thickness = zTn 186 pt%width = min(zWn,zLn)187 pt%length = max(zWn,zLn)189 pt%width = MIN( zWn , zLn ) 190 pt%length = MAX( zWn , zLn ) 188 191 189 192 next=>this%next … … 197 200 z1_e1e2 = r1_e1e2t(ii,ij) * this%mass_scaling 198 201 CALL icb_dia_size( ii, ij, zWn, zLn, zAbits, & 199 & this%mass_scaling, zMnew, znMbits, z1_e1e2)202 & this%mass_scaling, zMnew, znMbits, z1_e1e2 ) 200 203 ENDIF 201 204 ! … … 203 206 ! 204 207 END DO 205 208 206 209 ! now use melt and associated heat flux in ocean (or not) 207 210 !
Note: See TracChangeset
for help on using the changeset viewer.
