#include "cppdefs.h" MODULE ice_frazil_mod #if (defined ICE_MODEL && defined ICE_THERMO && !defined TS_FIXED) \ || defined CICE_MODEL ! !======================================================================= ! Copyright (c) 2002-2016 The ROMS/TOMS Group ! !================================================== Hernan G. Arango === ! ! ! This routine computes the frazil ice growth in the water when the ! water temperature gets below freezing. It adjusts the water ! temperature and salinity accordingly. ! ! Reference: Steele et al. (1989). JPO, 19, 139-147. ! ! !======================================================================= ! implicit none PRIVATE PUBLIC ice_frazil CONTAINS SUBROUTINE ice_frazil (ng, tile) USE mod_param USE mod_grid USE mod_ocean USE mod_ice USE mod_stepping # ifdef CICE_MODEL USE mod_forces # endif integer, intent(in) :: ng, tile ! # include "tile.h" ! # ifdef PROFILE CALL wclock_on (ng, iNLM, 51) # endif ! CALL ice_frazil_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & nnew(ng), & # ifdef MASKING & GRID(ng) % rmask, & # endif # ifdef WET_DRY & GRID(ng) % rmask_wet, & # endif & GRID(ng) % Hz, & & GRID(ng) % z_r, & & OCEAN(ng) % rho, & & OCEAN(ng) % t, & # ifdef CICE_MODEL & FORCES(ng) % sustr, & & FORCES(ng) % svstr, & & FORCES(ng) % stflx, & & ICE(ng) % ai, & & ICE(ng) % hi, & & GRID(ng) % z_w, & # endif & ICE(ng) % wfr) # ifdef PROFILE CALL wclock_off (ng, iNLM, 51) # endif RETURN END SUBROUTINE ice_frazil ! !*********************************************************************** subroutine ice_frazil_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & nnew, & # ifdef MASKING & rmask, & # endif # ifdef WET_DRY & rmask_wet, & # endif & Hz, z_r, rho, t, & # ifdef CICE_MODEL & sustr, svstr, stflx, ai, hi, z_w, & # endif & wfr) !*********************************************************************** ! USE mod_param USE mod_scalars # ifdef CICE_MODEL USE shr_const_mod # endif ! USE bc_2d_mod, ONLY : bc_r2d_tile USE exchange_3d_mod, ONLY : exchange_r3d_tile # ifdef DISTRIBUTE USE mp_exchange_mod USE distribute_mod, ONLY : mp_reduce # endif ! implicit none ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile integer, intent(in) :: LBi, UBi, LBj, UBj integer, intent(in) :: IminS, ImaxS, JminS, JmaxS integer, intent(in) :: nnew # ifdef ASSUMED_SHAPE # ifdef MASKING real(r8), intent(in) :: rmask(LBi:,LBj:) # endif # ifdef WET_DRY real(r8), intent(in) :: rmask_wet(LBi:,LBj:) # endif real(r8), intent(in) :: Hz(LBi:,LBj:,:) real(r8), intent(in) :: z_r(LBi:,LBj:,:) real(r8), intent(in) :: rho(LBi:,LBj:,:) real(r8), intent(inout) :: t(LBi:,LBj:,:,:,:) # ifdef CICE_MODEL real(r8), intent(in) :: sustr(LBi:,LBj:) real(r8), intent(in) :: svstr(LBi:,LBj:) real(r8), intent(in) :: stflx(LBi:,LBj:,:) real(r8), intent(in) :: ai(LBi:,LBj:) real(r8), intent(in) :: hi(LBi:,LBj:) real(r8), intent(in) :: z_w(LBi:,LBj:,0:) # endif real(r8), intent(out) :: wfr(LBi:,LBj:) # else # ifdef MASKING real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj) # endif # ifdef WET_DRY real(r8), intent(in) :: rmask_wet(LBi:UBi,LBj:UBj) # endif real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng)) real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng)) real(r8), intent(in) :: rho(LBi:UBi,LBj:UBj,N(ng)) real(r8), intent(inout) :: t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng)) # ifdef CICE_MODEL real(r8), intent(in) :: sustr(LBi:UBi,LBj:UBj) real(r8), intent(in) :: svstr(LBi:UBi,LBj:UBj) real(r8), intent(in) :: stflx(LBi:UBi,LBj:UBj,NT(ng)) real(r8), intent(in) :: ai(LBi:UBi,LBj:UBj) real(r8), intent(in) :: hi(LBi:UBi,LBj:UBj) real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng)) # endif real(r8), intent(out) :: wfr(LBi:UBi,LBj:UBj) # endif ! ! Local variable definitions ! integer :: i, j, k, itrc real(r8), parameter :: Lhat = 79.2_r8 real(r8), parameter :: r = 0.5_r8 real(r8) :: t_freeze real(r8) :: s1 real(r8) :: z1 real(r8) :: gamma_k real(r8) :: t_fr real(r8) :: ice_dens real(r8) :: delta_wfr # ifdef CICE_MODEL real(r8), parameter :: depressT = -0.054_r8 # else real(r8), parameter :: depressT = -0.0543_r8 # endif # ifdef CICE_MODEL # ifdef ICE_LOG_LAYER ! Compute heat flux from ocean a la Mellor and Kantha, 1989 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: t0mk real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: s0mk real(r8), parameter :: prt = 13._r8 real(r8), parameter :: prs = 2432._r8 real(r8), parameter :: tpr = 0.85_r8 real(r8), parameter :: nu = 1.8E-6_r8 real(r8), parameter :: z0ii = 0.02_r8 real(r8), parameter :: kappa = 0.4_r8 real(r8), parameter :: hfus = 3.347E+5_r8 ! [J kg-1] real(r8), parameter :: ykf = 3.14 real(r8), parameter :: eps = 1.0e-20_r8 real(r8) :: utau, dztop, rno, termt, terms, cht, chs, z0, zdz0 # endif real(r8) :: dz, heat, tfr real(r8), parameter :: bl_thick = 5.0_r8 real(r8) :: ice_thick, thickness, fract real(r8), parameter :: latent_heat_fusion = 3.337e5_r8 ! J/kg real(r8), parameter :: cp_sw = 3996.0_r8 ! J/kg/K real(r8), parameter :: rho_sw = 1026.0_r8 ! kg/m^3 real(r8), parameter :: rho_fw = 1000.0_r8 ! kg/m^3 real(r8), parameter :: cp_over_lhfusion = & & rho_sw*cp_sw/(latent_heat_fusion*rho_fw) real(r8), parameter :: rhocpr = 4093740._r8 ! [W s m-3 K-1] real(r8), parameter :: salice = 4.0_r8 real(r8), parameter :: salref = 34.7_r8 real(r8) :: potice # endif # ifdef DISTRIBUTE real(r8), allocatable :: buffer(:) character (len=3), allocatable :: op_handle(:) # endif ! Inline functions ! Freezing temperature (Gill, 1982) ! t_freeze(s1,z1) = -0.0575*s1 + 1.710523d-3*sqrt(s1)**3 ! & - 2.154996d-4*s1*s1 + 0.000753*z1 ! Freezing temperature (Steele et al. 1989) ! t_freeze(s1,z1) = -0.0543*s1 + 0.000759*z1 # ifdef CICE_MODEL ! Need to match CICE t_freeze(s1,z1) = -1.8_r8 ! t_freeze(s1,z1) = depressT*s1 # else ! Temperature or potential temperature?? t_freeze(s1,z1) = depressT*s1 # endif # include "set_bounds.h" # ifdef DISTRIBUTE IF (.not. allocated(buffer)) allocate (buffer(1)) IF (.not. allocated(op_handle)) allocate (op_handle(1)) # endif # ifdef CICE_MODEL ice_dens = SHR_CONST_RHOICE # else ice_dens = rhoice(ng) # endif DO j=Jstr,Jend DO i=Istr,Iend wfr(i,j) = 0.0_r8 ENDDO ENDDO # ifdef CICE_MODEL DO j=Jstr,Jend DO i=Istr,Iend # ifdef WET_DRY IF (rmask_wet(i,j) .ne. 0.0_r8) THEN # elif defined MASKING IF (rmask(i,j) .ne. 0.0_r8) THEN # endif thickness = 0.0_r8 DO k=N(ng),1,-1 dz = Hz(i,j,k) t_fr = t_freeze(t(i,j,k,nnew,isalt),z_r(i,j,k)) IF ((thickness + dz) > bl_thick) THEN fract = (bl_thick - thickness)/dz dz = fract*dz END IF thickness = thickness + dz potice = (t_fr - t(i,j,k,nnew,itemp)) * dz ! !*** if potice < 0, compute again later !*** if potice > 0, keep on freezing (wfr > 0) ! IF ((wfr(i,j) + potice) > 0) THEN wfr(i,j) = wfr(i,j) + potice ELSE potice = -wfr(i,j) wfr(i,j) = 0.0_r8 ENDIF t(i,j,k,nnew,itemp) = t(i,j,k,nnew,itemp) + & & potice/Hz(i,j,k) t(i,j,k,nnew,isalt) = t(i,j,k,nnew,isalt) & & + (salref-salice)*potice*cp_over_lhfusion/Hz(i,j,k) !*** jump out if done with surface layer IF (z_w(i,j,N(ng))-z_w(i,j,k) > bl_thick) EXIT END DO # if defined WET_DRY || defined MASKING END IF # endif END DO END DO # else ! Old ice model DO j=Jstr,Jend DO i=Istr,Iend DO k=1,N(ng) # ifdef WET_DRY IF (rmask_wet(i,j) .ne. 0.0_r8) THEN # elif defined MASKING IF (rmask(i,j) .ne. 0.0_r8) THEN # endif t_fr = t_freeze(t(i,j,k,nnew,isalt),z_r(i,j,k)) IF (t(i,j,k,nnew,itemp) .lt. t_fr) THEN gamma_k = (t_fr - t(i,j,k,nnew,itemp)) / & & (Lhat + t(i,j,k,nnew,itemp)*(1.0_r8 - r) & & - depressT * t(i,j,k,nnew,isalt)) IF (gamma_k .lt. 0.0_r8) THEN print *, 'trouble in ice_frazil', i, j, k, & & t(i,j,k,nnew,itemp), t(i,j,k,nnew,isalt), & & t_fr, wfr(i,j), gamma_k, Hz(i,j,k) exit_flag = 9 END IF wfr(i,j) = wfr(i,j) + gamma_k * Hz(i,j,k) * & & (rho0 + rho(i,j,k) ) / ice_dens t(i,j,k,nnew,itemp) = t(i,j,k,nnew,itemp) + gamma_k * & & (Lhat + t(i,j,k,nnew,itemp)*(1.0_r8 - r)) t(i,j,k,nnew,isalt) = t(i,j,k,nnew,isalt) * & & (1.0_r8 + gamma_k) ELSE IF (wfr(i,j) > 0 .and. & & t(i,j,k,nnew,itemp) .gt. t_fr) THEN ! Use heat at this level to melt some ice from below. ! gamma_k becomes negative here. gamma_k = (t_fr - t(i,j,k,nnew,itemp)) / & & (Lhat + t(i,j,k,nnew,itemp)*(1.0_r8 - r) & & - depressT * t(i,j,k,nnew,isalt)) delta_wfr = gamma_k * Hz(i,j,k) * & & (rho0 + rho(i,j,k) ) / rhoice(ng) IF ((wfr(i,j) + delta_wfr) > 0) THEN wfr(i,j) = wfr(i,j) + delta_wfr ELSE gamma_k = -wfr(i,j) * rhoice(ng) / & & (Hz(i,j,k)*(rho0+rho(i,j,k))) wfr(i,j) = 0.0_r8 ENDIF t(i,j,k,nnew,itemp) = t(i,j,k,nnew,itemp) + gamma_k * & & (Lhat + t(i,j,k,nnew,itemp)*(1.0_r8 - r)) t(i,j,k,nnew,isalt) = t(i,j,k,nnew,isalt) * & & (1.0_r8 + gamma_k) END IF # if defined WET_DRY || defined MASKING END IF # endif END DO wfr(i,j) = wfr(i,j)/dt(ng) IF (wfr(i,j) .lt. 0.0_r8) THEN print *, 'trouble in ice_frazil', i, j, & & t(i,j,N(ng),nnew,itemp), t(i,j,N(ng),nnew,isalt), & & wfr(i,j), gamma_k, Hz(i,j,N(ng)) exit_flag = 9 END IF END DO END DO # endif # ifdef CICE_MODEL ! Now compute heat flux from ocean and subtract from wfr. ! Using old salt flux, solve for surface salinity s0mk. # ifdef ICE_LOG_LAYER DO j=Jstr,Jend DO i=Istr,Iend dztop=z_w(i,j,N(ng))-z_r(i,j,N(ng)) ice_thick = 0.05_r8+hi(i,j)/MAX(ai(i,j),eps) utau = sqrt(sqrt( & & (0.5_r8*(sustr(i,j)+sustr(i+1,j)))**2 & & + (0.5_r8*(svstr(i,j)+svstr(i,j+1)))**2 & & ) ) utau = max(utau,1.E-4_r8) ! Need some roughness estimate here z0 = max(z0ii*ice_thick,0.01_r8) z0 = min(z0,0.1_r8) ! ! *** Yaglom and Kader formulation for z0t and z0s ! zdz0 = dztop/z0 !WPB zdz0 = MAX(zdz0,3._r8) ! Was this: rno = utau*0.09_r8/nu ! Mellor ! rno = utau*z0/nu termt = ykf*sqrt(rno)*prt**0.666667_r8 terms = ykf*sqrt(rno)*prs**0.666667_r8 cht = utau/(tpr*log(zdz0)/kappa+termt) chs = utau/(tpr*log(zdz0)/kappa+terms) IF (ai(i,j) .le. min_a(ng)) THEN s0mk(i,j) = t(i,j,N(ng),nnew,isalt) t0mk(i,j) = t(i,j,N(ng),nnew,itemp) ELSE s0mk(i,j) = t(i,j,N(ng),nnew,isalt) - & & stflx(i,j,isalt)/chs s0mk(i,j) = max(s0mk(i,j),0._r8) s0mk(i,j) = min(s0mk(i,j),40._r8) t0mk(i,j) = t_freeze(s0mk(i,j),0.0_r8) END IF ! convert units to W/m^2 wfr(i,j) = (wfr(i,j) & & -cht*(t(i,j,N(ng),nnew,itemp)-t0mk(i,j)))*rhocpr END DO END DO # else DO j=Jstr,Jend DO i=Istr,Iend heat = 0.0_r8 thickness = 0.0_r8 DO k=N(ng),1,-1 dz = Hz(i,j,k) tfr = t_freeze(t(i,j,k,nnew,isalt),0.0_r8) IF ((thickness + dz) < bl_thick) THEN thickness = thickness + dz heat = heat + (t(i,j,k,nnew,itemp)-tfr)*dz ELSE fract = (bl_thick - thickness)/dz ! thickness = thickness + dz*fract heat = heat + (t(i,j,k,nnew,itemp)-tfr)*dz*fract END IF IF (z_w(i,j,N(ng))-z_w(i,j,k) > bl_thick) EXIT END DO ! convert units to W/m^2 wfr(i,j) = (wfr(i,j) - heat)*rhocpr/dt(ng) END DO END DO # endif # endif CALL bc_r2d_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & wfr) # ifdef DISTRIBUTE buffer(1) = exit_flag op_handle(1) = 'MAX' CALL mp_reduce (ng, iNLM, 1, buffer, op_handle) exit_flag = int(buffer(1)) CALL mp_exchange2d (ng, tile, iNLM, 1, & & LBi, UBi, LBj, UBj, & & NghostPoints, EWperiodic(ng), NSperiodic(ng), & & wfr) # endif RETURN END SUBROUTINE ice_frazil_tile #endif END MODULE ice_frazil_mod