MODULE bdyice_lim !!====================================================================== !! *** MODULE bdyice_lim *** !! Unstructured Open Boundary Cond. : Open boundary conditions for sea-ice (LIM2 and LIM3) !!====================================================================== !! History : 3.3 ! 2010-09 (D. Storkey) Original code !! 3.4 ! 2011 (D. Storkey) rewrite in preparation for OBC-BDY merge !! - ! 2012-01 (C. Rousset) add lim3 and remove useless jk loop !!---------------------------------------------------------------------- #if defined key_bdy && ( defined key_lim2 || defined key_lim3 ) !!---------------------------------------------------------------------- !! 'key_bdy' and Unstructured Open Boundary Conditions !! 'key_lim2' LIM-2 sea ice model !! 'key_lim3' LIM-3 sea ice model !!---------------------------------------------------------------------- !! bdy_ice_lim : Application of open boundaries to ice !! bdy_ice_frs : Application of Flow Relaxation Scheme !!---------------------------------------------------------------------- USE timing ! Timing USE phycst ! physical constant USE eosbn2 ! equation of state USE oce ! ocean dynamics and tracers variables #if defined key_lim2 USE par_ice_2 USE ice_2 ! LIM_2 ice variables #elif defined key_lim3 USE par_ice USE ice ! LIM_3 ice variables #endif USE par_oce ! ocean parameters USE dom_oce ! ocean space and time domain variables USE dom_ice ! sea-ice domain USE sbc_oce ! Surface boundary condition: ocean fields USE bdy_oce ! ocean open boundary conditions USE lbclnk ! ocean lateral boundary conditions (or mpp link) USE in_out_manager ! write to numout file USE lib_mpp ! distributed memory computing USE lib_fortran ! to use key_nosignedzero IMPLICIT NONE PRIVATE PUBLIC bdy_ice_lim ! routine called in sbcmod PUBLIC bdy_ice_lim_dyn ! routine called in limrhg REAL(wp) :: epsi20 = 1.e-20_wp ! module constants !!---------------------------------------------------------------------- !! NEMO/OPA 3.3 , NEMO Consortium (2010) !! $Id: bdyice.F90 2715 2011-03-30 15:58:35Z rblod $ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE bdy_ice_lim( kt ) !!---------------------------------------------------------------------- !! *** SUBROUTINE bdy_ice_lim *** !! !! ** Purpose : - Apply open boundary conditions for ice (LIM2 and LIM3) !! !!---------------------------------------------------------------------- INTEGER, INTENT( in ) :: kt ! Main time step counter !! INTEGER :: ib_bdy ! Loop index DO ib_bdy=1, nb_bdy SELECT CASE( cn_ice_lim(ib_bdy) ) CASE('none') CYCLE CASE('frs') CALL bdy_ice_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy ) CASE DEFAULT CALL ctl_stop( 'bdy_ice_lim : unrecognised option for open boundaries for ice fields' ) END SELECT ENDDO END SUBROUTINE bdy_ice_lim SUBROUTINE bdy_ice_frs( idx, dta, kt, ib_bdy ) !!------------------------------------------------------------------------------ !! *** SUBROUTINE bdy_ice_frs *** !! !! ** Purpose : Apply the Flow Relaxation Scheme for sea-ice fields in the case !! of unstructured open boundaries. !! !! Reference : Engedahl H., 1995: Use of the flow relaxation scheme in a three- !! dimensional baroclinic ocean model with realistic topography. Tellus, 365-382. !!------------------------------------------------------------------------------ TYPE(OBC_INDEX), INTENT(in) :: idx ! OBC indices TYPE(OBC_DATA), INTENT(in) :: dta ! OBC external data INTEGER, INTENT(in) :: kt ! main time-step counter INTEGER, INTENT(in) :: ib_bdy ! BDY set index !! INTEGER :: jb, jk, jgrd, jl ! dummy loop indices INTEGER :: ji, jj ! local scalar REAL(wp) :: zwgt, zwgt1 ! local scalar REAL(wp) :: zinda, ztmelts !!------------------------------------------------------------------------------ ! IF( nn_timing == 1 ) CALL timing_start('bdy_ice_frs') ! !IF( kt /= nit000 ) THEN jgrd = 1 ! Everything is at T-points here ! #if defined key_lim2 DO jb = 1, idx%nblen(jgrd) ji = idx%nbi(jb,jgrd) jj = idx%nbj(jb,jgrd) zwgt = idx%nbw(jb,jgrd) zwgt1 = 1.e0 - idx%nbw(jb,jgrd) frld (ji,jj) = ( frld (ji,jj) * zwgt1 + dta%frld (jb) * zwgt ) * tmask(ji,jj,1) ! Leads fraction hicif(ji,jj) = ( hicif(ji,jj) * zwgt1 + dta%hicif(jb) * zwgt ) * tmask(ji,jj,1) ! Ice depth hsnif(ji,jj) = ( hsnif(ji,jj) * zwgt1 + dta%hsnif(jb) * zwgt ) * tmask(ji,jj,1) ! Snow depth ! zinda = 1.0 - MAX( 0.0_wp , SIGN ( 1.0_wp , - frld(ji,jj) ) ) ! 0 if no ice ! !------------------------------ ! ! Sea ice surface temperature ! !------------------------------ ! sist(ji,jj) = zinda * 270.0 + ( 1.0 - zinda ) * tfu(ji,jj) ! !----------------------------------------------- ! ! Ice/snow temperatures and energy stored in brines ! !----------------------------------------------- ! !!! TO BE CONTIUNED (as LIM3 below) !!! ! zindhe = MAX( 0.e0, SIGN( 1.e0, fcor(1,jj) ) ) ! = 0 for SH, =1 for NH ! ! ! Recover in situ values. ! zindb = MAX( rzero, SIGN( rone, zs0a(ji,jj) - epsi06 ) ) ! zacrith = 1.0 - ( zindhe * acrit(1) + ( 1.0 - zindhe ) * acrit(2) ) ! zs0a (ji,jj) = zindb * MIN( zs0a(ji,jj), zacrith ) ! hsnif(ji,jj) = zindb * ( zs0sn(ji,jj) /MAX( zs0a(ji,jj), epsi16 ) ) ! hicif(ji,jj) = zindb * ( zs0ice(ji,jj)/MAX( zs0a(ji,jj), epsi16 ) ) ! zindsn = MAX( rzero, SIGN( rone, hsnif(ji,jj) - epsi06 ) ) ! zindic = MAX( rzero, SIGN( rone, hicif(ji,jj) - epsi03 ) ) ! zindb = MAX( zindsn, zindic ) ! zs0a (ji,jj) = zindb * zs0a(ji,jj) ! frld (ji,jj) = 1.0 - zs0a(ji,jj) ! hsnif(ji,jj) = zindsn * hsnif(ji,jj) ! hicif(ji,jj) = zindic * hicif(ji,jj) ! zusvosn = 1.0/MAX( hsnif(ji,jj) * zs0a(ji,jj), epsi16 ) ! zusvoic = 1.0/MAX( hicif(ji,jj) * zs0a(ji,jj), epsi16 ) ! zignm = MAX( rzero, SIGN( rone, hsndif - hsnif(ji,jj) ) ) ! zrtt = 173.15 * rone ! ztsn = zignm * tbif(ji,jj,1) & ! + ( 1.0 - zignm ) * MIN( MAX( zrtt, rt0_snow * zusvosn * zs0c0(ji,jj)) , tfu(ji,jj) ) ! ztic1 = MIN( MAX( zrtt, rt0_ice * zusvoic * zs0c1(ji,jj) ) , tfu(ji,jj) ) ! ztic2 = MIN( MAX( zrtt, rt0_ice * zusvoic * zs0c2(ji,jj) ) , tfu(ji,jj) ) ! ! tbif(ji,jj,1) = zindsn * ztsn + ( 1.0 - zindsn ) * tfu(ji,jj) ! tbif(ji,jj,2) = zindic * ztic1 + ( 1.0 - zindic ) * tfu(ji,jj) ! tbif(ji,jj,3) = zindic * ztic2 + ( 1.0 - zindic ) * tfu(ji,jj) ! qstoif(ji,jj) = zindb * xlic * zs0st(ji,jj) / MAX( zs0a(ji,jj), epsi16 ) END DO CALL lbc_bdy_lnk( frld, 'T', 1., ib_bdy ) ! lateral boundary conditions CALL lbc_bdy_lnk( hicif, 'T', 1., ib_bdy ) CALL lbc_bdy_lnk( hsnif, 'T', 1., ib_bdy ) vt_i(:,:) = hicif(:,:) * frld(:,:) vt_s(:,:) = hsnif(:,:) * frld(:,:) ! #elif defined key_lim3 DO jl = 1, jpl DO jb = 1, idx%nblen(jgrd) ji = idx%nbi(jb,jgrd) jj = idx%nbj(jb,jgrd) zwgt = idx%nbw(jb,jgrd) zwgt1 = 1.e0 - idx%nbw(jb,jgrd) a_i (ji,jj,jl) = ( a_i (ji,jj,jl) * zwgt1 + dta%a_i (jb,jl) * zwgt ) * tmask(ji,jj,1) ! Leads fraction ht_i(ji,jj,jl) = ( ht_i(ji,jj,jl) * zwgt1 + dta%ht_i(jb,jl) * zwgt ) * tmask(ji,jj,1) ! Ice depth ht_s(ji,jj,jl) = ( ht_s(ji,jj,jl) * zwgt1 + dta%ht_s(jb,jl) * zwgt ) * tmask(ji,jj,1) ! Snow depth ENDDO CALL lbc_bdy_lnk( a_i(:,:,jl), 'T', 1., ib_bdy ) CALL lbc_bdy_lnk( ht_i(:,:,jl), 'T', 1., ib_bdy ) CALL lbc_bdy_lnk( ht_s(:,:,jl), 'T', 1., ib_bdy ) DO jb = 1, idx%nblen(jgrd) ji = idx%nbi(jb,jgrd) jj = idx%nbj(jb,jgrd) ! clem : condition on ice thickness depends on the ice velocity ! if velocity is outward (strictly), then ice thickness, volume... must be equal to adjacent values IF ( u_ice(ji+1,jj ) < 0. .AND. umask(ji-1,jj ,1) == 0. ) THEN ht_i(ji,jj,jl) = ht_i(ji+1,jj ,jl) a_i (ji,jj,jl) = a_i (ji+1,jj ,jl) ht_s(ji,jj,jl) = ht_s(ji+1,jj ,jl) ENDIF IF ( u_ice(ji-1,jj ) > 0. .AND. umask(ji+1,jj ,1) == 0. ) THEN ht_i(ji,jj,jl) = ht_i(ji-1,jj ,jl) a_i (ji,jj,jl) = a_i (ji-1,jj ,jl) ht_s(ji,jj,jl) = ht_s(ji-1,jj ,jl) ENDIF IF ( v_ice(ji ,jj+1) < 0. .AND. vmask(ji ,jj-1,1) == 0. ) THEN ht_i(ji,jj,jl) = ht_i(ji ,jj+1,jl) a_i (ji,jj,jl) = a_i (ji ,jj+1,jl) ht_s(ji,jj,jl) = ht_s(ji ,jj+1,jl) ENDIF IF ( v_ice(ji ,jj-1) > 0. .AND. vmask(ji ,jj+1,1) == 0. ) THEN ht_i(ji,jj,jl) = ht_i(ji ,jj-1,jl) a_i (ji,jj,jl) = a_i (ji ,jj-1,jl) ht_s(ji,jj,jl) = ht_s(ji ,jj-1,jl) ENDIF zinda = 1.0 - MAX( 0.0_wp , SIGN ( 1.0_wp , - a_i(ji,jj,jl) ) ) ! 0 if no ice ! -------------------- ! Ice and snow volumes ! -------------------- v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) v_s(ji,jj,jl) = ht_s(ji,jj,jl) * a_i(ji,jj,jl) ! ------------- ! Ice salinity !--------------- sm_i(ji,jj,jl) = zinda * bulk_sal + ( 1.0 - zinda ) * s_i_min smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) !---------- ! Ice age !---------- o_i(ji,jj,jl) = zinda * 1.0 + ( 1.0 - zinda ) oa_i(ji,jj,jl) = o_i(ji,jj,jl) * a_i(ji,jj,jl) !------------------------------ ! Sea ice surface temperature !------------------------------ t_su(ji,jj,jl) = zinda * 270.0 + ( 1.0 - zinda ) * t_bo(ji,jj) !------------------------------------ ! Snow temperature and heat content !------------------------------------ DO jk = 1, nlay_s t_s(ji,jj,jk,jl) = zinda * 270.00 + ( 1.0 - zinda ) * rtt ! Snow energy of melting e_s(ji,jj,jk,jl) = zinda * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus ) ! Change dimensions e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac ! Multiply by volume, so that heat content in 10^9 Joules e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * & v_s(ji,jj,jl) / nlay_s END DO !jk !----------------------------------------------- ! Ice salinities, temperature and heat content !----------------------------------------------- DO jk = 1, nlay_i t_i(ji,jj,jk,jl) = zinda * 270.00 + ( 1.0 - zinda ) * rtt s_i(ji,jj,jk,jl) = zinda * bulk_sal + ( 1.0 - zinda ) * s_i_min ztmelts = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K ! heat content per unit volume e_i(ji,jj,jk,jl) = zinda * rhoic * & ( cpic * ( ztmelts - t_i(ji,jj,jk,jl) ) & + lfus * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-epsi20) ) & - rcp * ( ztmelts - rtt ) ) ! Correct dimensions to avoid big values e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * & area(ji,jj) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / nlay_i END DO ! jk END DO !jb CALL lbc_bdy_lnk( a_i(:,:,jl), 'T', 1., ib_bdy ) ! lateral boundary conditions CALL lbc_bdy_lnk( ht_i(:,:,jl), 'T', 1., ib_bdy ) CALL lbc_bdy_lnk( ht_s(:,:,jl), 'T', 1., ib_bdy ) CALL lbc_bdy_lnk( v_i(:,:,jl), 'T', 1., ib_bdy ) CALL lbc_bdy_lnk( v_s(:,:,jl), 'T', 1., ib_bdy ) CALL lbc_bdy_lnk( smv_i(:,:,jl), 'T', 1., ib_bdy ) CALL lbc_bdy_lnk( sm_i(:,:,jl), 'T', 1., ib_bdy ) CALL lbc_bdy_lnk( oa_i(:,:,jl), 'T', 1., ib_bdy ) CALL lbc_bdy_lnk( o_i(:,:,jl), 'T', 1., ib_bdy ) CALL lbc_bdy_lnk( t_su(:,:,jl), 'T', 1., ib_bdy ) DO jk = 1, nlay_s CALL lbc_bdy_lnk(t_s(:,:,jk,jl), 'T', 1., ib_bdy ) CALL lbc_bdy_lnk(e_s(:,:,jk,jl), 'T', 1., ib_bdy ) END DO DO jk = 1, nlay_i CALL lbc_bdy_lnk(t_i(:,:,jk,jl), 'T', 1., ib_bdy ) CALL lbc_bdy_lnk(e_i(:,:,jk,jl), 'T', 1., ib_bdy ) END DO END DO !jl #endif ! !ENDIF !nit000/=0 IF( nn_timing == 1 ) CALL timing_stop('bdy_ice_frs') ! END SUBROUTINE bdy_ice_frs SUBROUTINE bdy_ice_lim_dyn( kn, pvar1, pvar2, pvar12 ) !!------------------------------------------------------------------------------ !! *** SUBROUTINE bdy_ice_lim_dyn *** !! !! ** Purpose : Apply dynamics boundary conditions for sea-ice in the cas of unstructured open boundaries. !! kn = 1: u_ice and v_ice are equal to the value of the adjacent grid point if this latter is not ice free !! if adjacent grid point is ice free, then u_ice and v_ice are equal to ocean velocities !! kn = 2: the stress tensor is set to 0 (i.e. pvar1, pvar2, pvar12) !! !! 2013-06 : C. Rousset !!------------------------------------------------------------------------------ !! INTEGER, INTENT( in ) :: kn ! set up of ice vel (kn=1) or stress tensor (kn=2) REAL(wp), INTENT( inout ), DIMENSION(:,:), OPTIONAL :: pvar1, pvar2, pvar12 ! stress tensor components (zs1,zs2,zs12) INTEGER :: jb, jgrd ! dummy loop indices INTEGER :: ji, jj ! local scalar INTEGER :: ib_bdy ! Loop index REAL(wp) :: zmsk1, zmsk2, zflag, zinda !!------------------------------------------------------------------------------ ! IF( nn_timing == 1 ) CALL timing_start('bdy_ice_lim_dyn') ! DO ib_bdy=1, nb_bdy ! SELECT CASE( nn_ice_lim(ib_bdy) ) CASE(jp_none) CYCLE CASE(jp_frs) IF ( kn == 1 ) THEN ! set up of u_ice and v_ice jgrd = 2 ! u velocity DO jb = 1, idx_bdy(ib_bdy)%nblen(jgrd) ji = idx_bdy(ib_bdy)%nbi(jb,jgrd) jj = idx_bdy(ib_bdy)%nbj(jb,jgrd) zflag = idx_bdy(ib_bdy)%flagu(jb) IF ( ABS( zflag ) == 1. ) THEN ! eastern and western boundaries ! one of the two zmsk is always 0 (because of zflag) zmsk1 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji+1,jj) ) ) ! 0 if no ice zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji-1,jj) ) ) ! 0 if no ice ! u_ice = u_ice of the adjacent grid point except if this grid point is ice-free (then u_ice = u_oce) u_ice (ji,jj) = u_ice(ji+1,jj) * 0.5 * ABS( zflag + 1._wp ) * zmsk1 + & & u_ice(ji-1,jj) * 0.5 * ABS( zflag - 1._wp ) * zmsk2 + & & u_oce(ji ,jj) * ( 1._wp - MIN( 1._wp, zmsk1 + zmsk2 ) ) ELSE ! everywhere else u_ice(ji,jj) = u_oce(ji,jj) ENDIF ! mask ice velocities zinda = 1.0 - MAX( 0.0_wp, SIGN( 1.0_wp , - at_i(ji,jj) ) ) ! 0 if no ice u_ice(ji,jj) = zinda * u_ice(ji,jj) ENDDO jgrd = 3 ! v velocity DO jb = 1, idx_bdy(ib_bdy)%nblen(jgrd) ji = idx_bdy(ib_bdy)%nbi(jb,jgrd) jj = idx_bdy(ib_bdy)%nbj(jb,jgrd) zflag = idx_bdy(ib_bdy)%flagv(jb) IF ( ABS( zflag ) == 1. ) THEN ! northern and southern boundaries ! one of the two zmsk is always 0 (because of zflag) zmsk1 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj+1) ) ) ! 0 if no ice zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj-1) ) ) ! 0 if no ice ! u_ice = u_ice of the adjacent grid point except if this grid point is ice-free (then u_ice = u_oce) v_ice (ji,jj) = v_ice(ji,jj+1) * 0.5 * ABS( zflag + 1._wp ) * zmsk1 + & & v_ice(ji,jj-1) * 0.5 * ABS( zflag - 1._wp ) * zmsk2 + & & v_oce(ji,jj ) * ( 1._wp - MIN( 1._wp, zmsk1 + zmsk2 ) ) ELSE ! everywhere else v_ice(ji,jj) = v_oce(ji,jj) ENDIF ! mask ice velocities zinda = 1.0 - MAX( 0.0_wp, SIGN( 1.0_wp , - at_i(ji,jj) ) ) ! 0 if no ice v_ice(ji,jj) = zinda * v_ice(ji,jj) ENDDO CALL lbc_bdy_lnk( u_ice(:,:), 'U', -1., ib_bdy ) CALL lbc_bdy_lnk( v_ice(:,:), 'V', -1., ib_bdy ) ELSEIF ( kn == 2 ) THEN ! set up of stress tensor (not sure it works) jgrd = 1 ! T grid DO jb = 1, idx_bdy(ib_bdy)%nblen(jgrd) ji = idx_bdy(ib_bdy)%nbi(jb,jgrd) jj = idx_bdy(ib_bdy)%nbj(jb,jgrd) pvar1 (ji,jj) = 0._wp pvar2 (ji,jj) = 0._wp pvar12(ji,jj) = 0._wp ENDDO ENDIF CASE DEFAULT CALL ctl_stop( 'bdy_ice_lim_dyn : unrecognised option for open boundaries for ice fields' ) END SELECT ENDDO IF( nn_timing == 1 ) CALL timing_stop('bdy_ice_lim_dyn') END SUBROUTINE bdy_ice_lim_dyn #else !!--------------------------------------------------------------------------------- !! Default option Empty module !!--------------------------------------------------------------------------------- CONTAINS SUBROUTINE bdy_ice_lim( kt ) ! Empty routine WRITE(*,*) 'bdy_ice_lim: You should not have seen this print! error?', kt END SUBROUTINE bdy_ice_lim #endif !!================================================================================= END MODULE bdyice_lim