- Timestamp:
- 2018-11-07T18:25:49+01:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2018/dev_r9866_HPC_03_globcom/src/OCE/BDY/bdyice.F90
- Property svn:keywords set to Id
r9657 r10288 18 18 USE ice ! sea-ice: variables 19 19 USE icevar ! sea-ice: operations 20 USE ice itd ! sea-ice: rebining20 USE icecor ! sea-ice: corrections 21 21 USE icectl ! sea-ice: control prints 22 22 USE phycst ! physical constant … … 41 41 !!---------------------------------------------------------------------- 42 42 !! NEMO/OCE 4.0 , NEMO Consortium (2018) 43 !! $Id : bdyice.F90 8306 2017-07-10 10:18:03Z clem$44 !! Software governed by the CeCILL licen ce (./LICENSE)43 !! $Id$ 44 !! Software governed by the CeCILL license (see ./LICENSE) 45 45 !!---------------------------------------------------------------------- 46 46 CONTAINS … … 50 50 !! *** SUBROUTINE bdy_ice *** 51 51 !! 52 !! ** Purpose : - Apply open boundary conditions for ice (SI3)52 !! ** Purpose : Apply open boundary conditions for sea ice 53 53 !! 54 54 !!---------------------------------------------------------------------- 55 55 INTEGER, INTENT(in) :: kt ! Main time step counter 56 56 ! 57 INTEGER :: ib_bdy ! Loopindex57 INTEGER :: jbdy ! BDY set index 58 58 !!---------------------------------------------------------------------- 59 59 ! 60 IF( ln_timing ) CALL timing_start('bdy_ice ')60 IF( ln_timing ) CALL timing_start('bdy_ice_thd') 61 61 ! 62 62 CALL ice_var_glo2eqv 63 63 ! 64 DO ib_bdy = 1, nb_bdy65 ! 66 SELECT CASE( cn_ice( ib_bdy) )64 DO jbdy = 1, nb_bdy 65 ! 66 SELECT CASE( cn_ice(jbdy) ) 67 67 CASE('none') ; CYCLE 68 CASE('frs' ) ; CALL bdy_ice_frs( idx_bdy( ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )68 CASE('frs' ) ; CALL bdy_ice_frs( idx_bdy(jbdy), dta_bdy(jbdy), kt, jbdy ) 69 69 CASE DEFAULT 70 70 CALL ctl_stop( 'bdy_ice : unrecognised option for open boundaries for ice fields' ) … … 73 73 END DO 74 74 ! 75 CALL ice_var_zapsmall 75 CALL ice_cor( kt , 0 ) ! -- In case categories are out of bounds, do a remapping 76 ! ! i.e. inputs have not the same ice thickness distribution (set by rn_himean) 77 ! ! than the regional simulation 76 78 CALL ice_var_agg(1) 77 79 ! 78 80 IF( ln_icectl ) CALL ice_prt( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' ) 79 IF( ln_timing ) CALL timing_stop('bdy_ice ')81 IF( ln_timing ) CALL timing_stop('bdy_ice_thd') 80 82 ! 81 83 END SUBROUTINE bdy_ice 82 84 83 85 84 SUBROUTINE bdy_ice_frs( idx, dta, kt, ib_bdy )86 SUBROUTINE bdy_ice_frs( idx, dta, kt, jbdy ) 85 87 !!------------------------------------------------------------------------------ 86 88 !! *** SUBROUTINE bdy_ice_frs *** 87 89 !! 88 !! ** Purpose : Apply the Flow Relaxation Scheme for sea-ice fields in the case 89 !! of unstructured open boundaries. 90 !! ** Purpose : Apply the Flow Relaxation Scheme for sea-ice fields 90 91 !! 91 92 !! Reference : Engedahl H., 1995: Use of the flow relaxation scheme in a three- … … 95 96 TYPE(OBC_DATA), INTENT(in) :: dta ! OBC external data 96 97 INTEGER, INTENT(in) :: kt ! main time-step counter 97 INTEGER, INTENT(in) :: ib_bdy! BDY set index98 INTEGER, INTENT(in) :: jbdy ! BDY set index 98 99 ! 99 100 INTEGER :: jpbound ! 0 = incoming ice 100 101 ! ! 1 = outgoing ice 101 INTEGER :: jb, jk, jgrd, jl! dummy loop indices102 INTEGER :: ji, jj, ii, ij ! local scalar102 INTEGER :: i_bdy, jgrd ! dummy loop indices 103 INTEGER :: ji, jj, jk, jl, ib, jb 103 104 REAL(wp) :: zwgt, zwgt1 ! local scalar 104 105 REAL(wp) :: ztmelts, zdh … … 108 109 ! 109 110 DO jl = 1, jpl 110 DO jb= 1, idx%nblenrim(jgrd)111 ji = idx%nbi( jb,jgrd)112 jj = idx%nbj( jb,jgrd)113 zwgt = idx%nbw( jb,jgrd)114 zwgt1 = 1.e0 - idx%nbw( jb,jgrd)115 a_i(ji,jj,jl) = ( a_i(ji,jj,jl) * zwgt1 + dta%a_i( jb,jl) * zwgt ) * tmask(ji,jj,1) ! Leads fraction116 h_i(ji,jj,jl) = ( h_i(ji,jj,jl) * zwgt1 + dta%h_i( jb,jl) * zwgt ) * tmask(ji,jj,1) ! Ice depth117 h_s(ji,jj,jl) = ( h_s(ji,jj,jl) * zwgt1 + dta%h_s( jb,jl) * zwgt ) * tmask(ji,jj,1) ! Snow depth111 DO i_bdy = 1, idx%nblenrim(jgrd) 112 ji = idx%nbi(i_bdy,jgrd) 113 jj = idx%nbj(i_bdy,jgrd) 114 zwgt = idx%nbw(i_bdy,jgrd) 115 zwgt1 = 1.e0 - idx%nbw(i_bdy,jgrd) 116 a_i(ji,jj,jl) = ( a_i(ji,jj,jl) * zwgt1 + dta%a_i(i_bdy,jl) * zwgt ) * tmask(ji,jj,1) ! Leads fraction 117 h_i(ji,jj,jl) = ( h_i(ji,jj,jl) * zwgt1 + dta%h_i(i_bdy,jl) * zwgt ) * tmask(ji,jj,1) ! Ice depth 118 h_s(ji,jj,jl) = ( h_s(ji,jj,jl) * zwgt1 + dta%h_s(i_bdy,jl) * zwgt ) * tmask(ji,jj,1) ! Snow depth 118 119 119 120 ! ----------------- … … 124 125 125 126 ! Then, a) transfer the snow excess into the ice (different from icethd_dh) 126 zdh = MAX( 0._wp, ( rhos n * h_s(ji,jj,jl) + ( rhoic- rau0 ) * h_i(ji,jj,jl) ) * r1_rau0 )127 zdh = MAX( 0._wp, ( rhos * h_s(ji,jj,jl) + ( rhoi - rau0 ) * h_i(ji,jj,jl) ) * r1_rau0 ) 127 128 ! Or, b) transfer all the snow into ice (if incoming ice is likely to melt as it comes into a warmer environment) 128 !zdh = MAX( 0._wp, h_s(ji,jj,jl) * rhos n / rhoic)129 !zdh = MAX( 0._wp, h_s(ji,jj,jl) * rhos / rhoi ) 129 130 130 131 ! recompute h_i, h_s 131 132 h_i(ji,jj,jl) = MIN( hi_max(jl), h_i(ji,jj,jl) + zdh ) 132 h_s(ji,jj,jl) = MAX( 0._wp, h_s(ji,jj,jl) - zdh * rhoi c / rhosn)133 h_s(ji,jj,jl) = MAX( 0._wp, h_s(ji,jj,jl) - zdh * rhoi / rhos ) 133 134 134 135 ENDDO 135 CALL lbc_bdy_lnk( a_i(:,:,jl), 'T', 1., ib_bdy )136 CALL lbc_bdy_lnk( h_i(:,:,jl), 'T', 1., ib_bdy )137 CALL lbc_bdy_lnk( h_s(:,:,jl), 'T', 1., ib_bdy )138 136 ENDDO 139 ! retrieve at_i 140 at_i(:,:) = 0._wp 137 CALL lbc_bdy_lnk( a_i(:,:,:), 'T', 1., jbdy ) 138 CALL lbc_bdy_lnk( h_i(:,:,:), 'T', 1., jbdy ) 139 CALL lbc_bdy_lnk( h_s(:,:,:), 'T', 1., jbdy ) 140 141 141 DO jl = 1, jpl 142 at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 143 END DO 144 145 DO jl = 1, jpl 146 DO jb = 1, idx%nblenrim(jgrd) 147 ji = idx%nbi(jb,jgrd) 148 jj = idx%nbj(jb,jgrd) 142 DO i_bdy = 1, idx%nblenrim(jgrd) 143 ji = idx%nbi(i_bdy,jgrd) 144 jj = idx%nbj(i_bdy,jgrd) 149 145 150 146 ! condition on ice thickness depends on the ice velocity 151 147 ! if velocity is outward (strictly), then ice thickness, volume... must be equal to adjacent values 152 jpbound = 0 ; ii = ji ; ij = jj 153 ! 154 IF( u_ice(ji+1,jj ) < 0. .AND. umask(ji-1,jj ,1) == 0. ) jpbound = 1; ii = ji+1; ij = jj 155 IF( u_ice(ji-1,jj ) > 0. .AND. umask(ji+1,jj ,1) == 0. ) jpbound = 1; ii = ji-1; ij = jj 156 IF( v_ice(ji ,jj+1) < 0. .AND. vmask(ji ,jj-1,1) == 0. ) jpbound = 1; ii = ji ; ij = jj+1 157 IF( v_ice(ji ,jj-1) > 0. .AND. vmask(ji ,jj+1,1) == 0. ) jpbound = 1; ii = ji ; ij = jj-1 158 ! 159 IF( nn_ice_dta(ib_bdy) == 0 ) jpbound = 0; ii = ji; ij = jj ! case ice boundaries = initial conditions 160 ! ! do not make state variables dependent on velocity 161 ! 162 rswitch = MAX( 0.0_wp , SIGN ( 1.0_wp , at_i(ii,ij) - 0.01 ) ) ! 0 if no ice 163 ! 164 ! concentration and thickness 165 a_i(ji,jj,jl) = a_i(ii,ij,jl) * rswitch 166 h_i(ji,jj,jl) = h_i(ii,ij,jl) * rswitch 167 h_s(ji,jj,jl) = h_s(ii,ij,jl) * rswitch 168 ! 169 ! Ice and snow volumes 170 v_i(ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl) 171 v_s(ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl) 172 ! 173 SELECT CASE( jpbound ) 174 ! 175 CASE( 0 ) ! velocity is inward 176 ! 177 ! Ice salinity, age, temperature 178 s_i (ji,jj,jl) = rswitch * rn_ice_sal(ib_bdy) + ( 1.0 - rswitch ) * rn_simin 179 oa_i(ji,jj,jl) = rswitch * rn_ice_age(ib_bdy) * a_i(ji,jj,jl) 180 t_su(ji,jj,jl) = rswitch * rn_ice_tem(ib_bdy) + ( 1.0 - rswitch ) * rn_ice_tem(ib_bdy) 148 jpbound = 0 ; ib = ji ; jb = jj 149 ! 150 IF( u_ice(ji+1,jj ) < 0. .AND. umask(ji-1,jj ,1) == 0. ) jpbound = 1 ; ib = ji+1 ; jb = jj 151 IF( u_ice(ji-1,jj ) > 0. .AND. umask(ji+1,jj ,1) == 0. ) jpbound = 1 ; ib = ji-1 ; jb = jj 152 IF( v_ice(ji ,jj+1) < 0. .AND. vmask(ji ,jj-1,1) == 0. ) jpbound = 1 ; ib = ji ; jb = jj+1 153 IF( v_ice(ji ,jj-1) > 0. .AND. vmask(ji ,jj+1,1) == 0. ) jpbound = 1 ; ib = ji ; jb = jj-1 154 ! 155 IF( nn_ice_dta(jbdy) == 0 ) jpbound = 0 ; ib = ji ; jb = jj ! case ice boundaries = initial conditions 156 ! ! do not make state variables dependent on velocity 157 ! 158 IF( a_i(ib,jb,jl) > 0._wp ) THEN ! there is ice at the boundary 159 ! 160 a_i(ji,jj,jl) = a_i(ib,jb,jl) ! concentration 161 h_i(ji,jj,jl) = h_i(ib,jb,jl) ! thickness ice 162 h_s(ji,jj,jl) = h_s(ib,jb,jl) ! thickness snw 163 ! 164 SELECT CASE( jpbound ) 165 ! 166 CASE( 0 ) ! velocity is inward 167 ! 168 oa_i(ji,jj, jl) = rn_ice_age(jbdy) * a_i(ji,jj,jl) ! age 169 a_ip(ji,jj, jl) = 0._wp ! pond concentration 170 v_ip(ji,jj, jl) = 0._wp ! pond volume 171 t_su(ji,jj, jl) = rn_ice_tem(jbdy) ! temperature surface 172 t_s (ji,jj,:,jl) = rn_ice_tem(jbdy) ! temperature snw 173 t_i (ji,jj,:,jl) = rn_ice_tem(jbdy) ! temperature ice 174 s_i (ji,jj, jl) = rn_ice_sal(jbdy) ! salinity 175 sz_i(ji,jj,:,jl) = rn_ice_sal(jbdy) ! salinity profile 176 ! 177 CASE( 1 ) ! velocity is outward 178 ! 179 oa_i(ji,jj, jl) = oa_i(ib,jb, jl) ! age 180 a_ip(ji,jj, jl) = a_ip(ib,jb, jl) ! pond concentration 181 v_ip(ji,jj, jl) = v_ip(ib,jb, jl) ! pond volume 182 t_su(ji,jj, jl) = t_su(ib,jb, jl) ! temperature surface 183 t_s (ji,jj,:,jl) = t_s (ib,jb,:,jl) ! temperature snw 184 t_i (ji,jj,:,jl) = t_i (ib,jb,:,jl) ! temperature ice 185 s_i (ji,jj, jl) = s_i (ib,jb, jl) ! salinity 186 sz_i(ji,jj,:,jl) = sz_i(ib,jb,:,jl) ! salinity profile 187 ! 188 END SELECT 189 ! 190 IF( nn_icesal == 1 ) THEN ! if constant salinity 191 s_i (ji,jj ,jl) = rn_icesal 192 sz_i(ji,jj,:,jl) = rn_icesal 193 ENDIF 194 ! 195 ! global fields 196 v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl) ! volume ice 197 v_s (ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl) ! volume snw 198 sv_i(ji,jj,jl) = MIN( s_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content 181 199 DO jk = 1, nlay_s 182 t_s(ji,jj,jk,jl) = rswitch * rn_ice_tem(ib_bdy) + ( 1.0 - rswitch ) * rt0 200 e_s(ji,jj,jk,jl) = rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus ) ! enthalpy in J/m3 201 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s ! enthalpy in J/m2 202 END DO 203 DO jk = 1, nlay_i 204 ztmelts = - rTmlt * sz_i(ji,jj,jk,jl) ! Melting temperature in C 205 t_i(ji,jj,jk,jl) = MIN( t_i(ji,jj,jk,jl), ztmelts + rt0 ) ! Force t_i to be lower than melting point => likely conservation issue 206 ! 207 e_i(ji,jj,jk,jl) = rhoi * ( rcpi * ( ztmelts - ( t_i(ji,jj,jk,jl) - rt0 ) ) & ! enthalpy in J/m3 208 & + rLfus * ( 1._wp - ztmelts / ( t_i(ji,jj,jk,jl) - rt0 ) ) & 209 & - rcp * ztmelts ) 210 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i ! enthalpy in J/m2 183 211 END DO 184 DO jk = 1, nlay_i 185 t_i (ji,jj,jk,jl) = rswitch * rn_ice_tem(ib_bdy) + ( 1.0 - rswitch ) * rt0 186 sz_i(ji,jj,jk,jl) = rswitch * rn_ice_sal(ib_bdy) + ( 1.0 - rswitch ) * rn_simin 187 END DO 188 ! 189 CASE( 1 ) ! velocity is outward 190 ! 191 ! Ice salinity, age, temperature 192 s_i (ji,jj,jl) = rswitch * s_i (ii,ij,jl) + ( 1.0 - rswitch ) * rn_simin 193 oa_i(ji,jj,jl) = rswitch * oa_i(ii,ij,jl) 194 t_su(ji,jj,jl) = rswitch * t_su(ii,ij,jl) + ( 1.0 - rswitch ) * rt0 195 DO jk = 1, nlay_s 196 t_s(ji,jj,jk,jl) = rswitch * t_s(ii,ij,jk,jl) + ( 1.0 - rswitch ) * rt0 197 END DO 198 DO jk = 1, nlay_i 199 t_i (ji,jj,jk,jl) = rswitch * t_i (ii,ij,jk,jl) + ( 1.0 - rswitch ) * rt0 200 sz_i(ji,jj,jk,jl) = rswitch * sz_i(ii,ij,jk,jl) + ( 1.0 - rswitch ) * rn_simin 201 END DO 202 ! 203 END SELECT 204 ! 205 IF( nn_icesal == 1 ) THEN ! constant salinity : overwrite rn_icesal 206 s_i (ji,jj ,jl) = rn_icesal 207 sz_i(ji,jj,:,jl) = rn_icesal 212 ! 213 ELSE ! no ice at the boundary 214 ! 215 a_i (ji,jj, jl) = 0._wp 216 h_i (ji,jj, jl) = 0._wp 217 h_s (ji,jj, jl) = 0._wp 218 oa_i(ji,jj, jl) = 0._wp 219 a_ip(ji,jj, jl) = 0._wp 220 v_ip(ji,jj, jl) = 0._wp 221 t_su(ji,jj, jl) = rt0 222 t_s (ji,jj,:,jl) = rt0 223 t_i (ji,jj,:,jl) = rt0 224 225 IF( nn_icesal == 1 ) THEN ! if constant salinity 226 s_i (ji,jj ,jl) = rn_icesal 227 sz_i(ji,jj,:,jl) = rn_icesal 228 ELSE ! if variable salinity 229 s_i (ji,jj,jl) = rn_simin 230 sz_i(ji,jj,:,jl) = rn_simin 231 ENDIF 232 ! 233 ! global fields 234 v_i (ji,jj, jl) = 0._wp 235 v_s (ji,jj, jl) = 0._wp 236 sv_i(ji,jj, jl) = 0._wp 237 e_s (ji,jj,:,jl) = 0._wp 238 e_i (ji,jj,:,jl) = 0._wp 239 208 240 ENDIF 209 ! 210 ! contents 211 sv_i(ji,jj,jl) = MIN( s_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) 212 DO jk = 1, nlay_s 213 ! Snow energy of melting 214 e_s(ji,jj,jk,jl) = rswitch * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus ) 215 ! Multiply by volume, so that heat content in J/m2 216 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s 217 END DO 218 DO jk = 1, nlay_i 219 ztmelts = - tmut * sz_i(ji,jj,jk,jl) + rt0 !Melting temperature in K 220 ! heat content per unit volume 221 e_i(ji,jj,jk,jl) = rswitch * rhoic * & 222 ( cpic * ( ztmelts - t_i(ji,jj,jk,jl) ) & 223 + lfus * ( 1.0 - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) & 224 - rcp * ( ztmelts - rt0 ) ) 225 ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2 226 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * a_i(ji,jj,jl) * h_i(ji,jj,jl) * r1_nlay_i 227 END DO 228 ! 241 229 242 END DO 230 243 ! 231 CALL lbc_bdy_lnk( a_i(:,:,jl), 'T', 1., ib_bdy ) 232 CALL lbc_bdy_lnk( h_i(:,:,jl), 'T', 1., ib_bdy ) 233 CALL lbc_bdy_lnk( h_s(:,:,jl), 'T', 1., ib_bdy ) 234 CALL lbc_bdy_lnk( v_i(:,:,jl), 'T', 1., ib_bdy ) 235 CALL lbc_bdy_lnk( v_s(:,:,jl), 'T', 1., ib_bdy ) 236 ! 237 CALL lbc_bdy_lnk( sv_i(:,:,jl), 'T', 1., ib_bdy ) 238 CALL lbc_bdy_lnk( s_i(:,:,jl), 'T', 1., ib_bdy ) 239 CALL lbc_bdy_lnk( oa_i(:,:,jl), 'T', 1., ib_bdy ) 240 CALL lbc_bdy_lnk( t_su(:,:,jl), 'T', 1., ib_bdy ) 241 DO jk = 1, nlay_s 242 CALL lbc_bdy_lnk(t_s(:,:,jk,jl), 'T', 1., ib_bdy ) 243 CALL lbc_bdy_lnk(e_s(:,:,jk,jl), 'T', 1., ib_bdy ) 244 END DO 245 DO jk = 1, nlay_i 246 CALL lbc_bdy_lnk(t_i(:,:,jk,jl), 'T', 1., ib_bdy ) 247 CALL lbc_bdy_lnk(e_i(:,:,jk,jl), 'T', 1., ib_bdy ) 248 END DO 249 ! 250 END DO !jl 251 ! 252 ! --- In case categories are out of bounds, do a remapping --- ! 253 ! i.e. inputs have not the same ice thickness distribution 254 ! (set by rn_himean) than the regional simulation 255 IF( jpl > 1 ) CALL ice_itd_reb( kt ) 244 END DO ! jl 245 246 CALL lbc_bdy_lnk( a_i (:,:,:) , 'T', 1., jbdy ) 247 CALL lbc_bdy_lnk( h_i (:,:,:) , 'T', 1., jbdy ) 248 CALL lbc_bdy_lnk( h_s (:,:,:) , 'T', 1., jbdy ) 249 CALL lbc_bdy_lnk( oa_i(:,:,:) , 'T', 1., jbdy ) 250 CALL lbc_bdy_lnk( a_ip(:,:,:) , 'T', 1., jbdy ) 251 CALL lbc_bdy_lnk( v_ip(:,:,:) , 'T', 1., jbdy ) 252 CALL lbc_bdy_lnk( s_i (:,:,:) , 'T', 1., jbdy ) 253 CALL lbc_bdy_lnk( t_su(:,:,:) , 'T', 1., jbdy ) 254 CALL lbc_bdy_lnk( v_i (:,:,:) , 'T', 1., jbdy ) 255 CALL lbc_bdy_lnk( v_s (:,:,:) , 'T', 1., jbdy ) 256 CALL lbc_bdy_lnk( sv_i(:,:,:) , 'T', 1., jbdy ) 257 CALL lbc_bdy_lnk( t_s (:,:,:,:), 'T', 1., jbdy ) 258 CALL lbc_bdy_lnk( e_s (:,:,:,:), 'T', 1., jbdy ) 259 CALL lbc_bdy_lnk( t_i (:,:,:,:), 'T', 1., jbdy ) 260 CALL lbc_bdy_lnk( e_i (:,:,:,:), 'T', 1., jbdy ) 256 261 ! 257 262 END SUBROUTINE bdy_ice_frs … … 262 267 !! *** SUBROUTINE bdy_ice_dyn *** 263 268 !! 264 !! ** Purpose : Apply dynamics boundary conditions for sea-ice in the cas of unstructured open boundaries. 265 !! u_ice and v_ice are equal to the value of the adjacent grid point if this latter is not ice free 266 !! if adjacent grid point is ice free, then u_ice and v_ice are equal to ocean velocities 269 !! ** Purpose : Apply dynamics boundary conditions for sea-ice. 267 270 !! 268 !! 2013-06 : C. Rousset 271 !! ** Method : if this adjacent grid point is not ice free, then u_ice and v_ice take its value 272 !! if is ice free, then u_ice and v_ice are unchanged by BDY 273 !! they keep values calculated in rheology 274 !! 269 275 !!------------------------------------------------------------------------------ 270 276 CHARACTER(len=1), INTENT(in) :: cd_type ! nature of velocity grid-points 271 277 ! 272 INTEGER :: jb, jgrd! dummy loop indices273 INTEGER :: ji, jj 274 INTEGER :: ib_bdy ! Loopindex278 INTEGER :: i_bdy, jgrd ! dummy loop indices 279 INTEGER :: ji, jj ! local scalar 280 INTEGER :: jbdy ! BDY set index 275 281 REAL(wp) :: zmsk1, zmsk2, zflag 276 282 !!------------------------------------------------------------------------------ 277 ! 278 DO ib_bdy=1, nb_bdy 279 ! 280 SELECT CASE( cn_ice(ib_bdy) ) 283 IF( ln_timing ) CALL timing_start('bdy_ice_dyn') 284 ! 285 DO jbdy=1, nb_bdy 286 ! 287 SELECT CASE( cn_ice(jbdy) ) 281 288 ! 282 289 CASE('none') … … 285 292 CASE('frs') 286 293 ! 287 IF( nn_ice_dta( ib_bdy) == 0 ) CYCLE ! case ice boundaries = initial conditions288 ! 294 IF( nn_ice_dta(jbdy) == 0 ) CYCLE ! case ice boundaries = initial conditions 295 ! ! do not change ice velocity (it is only computed by rheology) 289 296 SELECT CASE ( cd_type ) 290 297 ! 291 298 CASE ( 'U' ) 292 299 jgrd = 2 ! u velocity 293 DO jb = 1, idx_bdy(ib_bdy)%nblenrim(jgrd)294 ji = idx_bdy( ib_bdy)%nbi(jb,jgrd)295 jj = idx_bdy( ib_bdy)%nbj(jb,jgrd)296 zflag = idx_bdy( ib_bdy)%flagu(jb,jgrd)300 DO i_bdy = 1, idx_bdy(jbdy)%nblenrim(jgrd) 301 ji = idx_bdy(jbdy)%nbi(i_bdy,jgrd) 302 jj = idx_bdy(jbdy)%nbj(i_bdy,jgrd) 303 zflag = idx_bdy(jbdy)%flagu(i_bdy,jgrd) 297 304 ! 298 305 IF ( ABS( zflag ) == 1. ) THEN ! eastern and western boundaries … … 301 308 zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji-1,jj) ) ) ! 0 if no ice 302 309 ! 303 ! u_ice = u_ice of the adjacent grid point except if this grid point is ice-free (then u_ice = u_oce)310 ! u_ice = u_ice of the adjacent grid point except if this grid point is ice-free (then do not change u_ice) 304 311 u_ice (ji,jj) = u_ice(ji+1,jj) * 0.5_wp * ABS( zflag + 1._wp ) * zmsk1 + & 305 312 & u_ice(ji-1,jj) * 0.5_wp * ABS( zflag - 1._wp ) * zmsk2 + & 306 & u_ oce(ji ,jj) * ( 1._wp - MIN( 1._wp, zmsk1 + zmsk2 ) )313 & u_ice(ji ,jj) * ( 1._wp - MIN( 1._wp, zmsk1 + zmsk2 ) ) 307 314 ELSE ! everywhere else 308 !u_ice(ji,jj) = u_oce(ji,jj)309 315 u_ice(ji,jj) = 0._wp 310 316 ENDIF 311 ! mask ice velocities312 rswitch = MAX( 0.0_wp , SIGN ( 1.0_wp , at_i(ji,jj) - 0.01_wp ) ) ! 0 if no ice313 u_ice(ji,jj) = rswitch * u_ice(ji,jj)314 317 ! 315 318 END DO 316 CALL lbc_bdy_lnk( u_ice(:,:), 'U', -1., ib_bdy )319 CALL lbc_bdy_lnk( u_ice(:,:), 'U', -1., jbdy ) 317 320 ! 318 321 CASE ( 'V' ) 319 322 jgrd = 3 ! v velocity 320 DO jb = 1, idx_bdy(ib_bdy)%nblenrim(jgrd)321 ji = idx_bdy( ib_bdy)%nbi(jb,jgrd)322 jj = idx_bdy( ib_bdy)%nbj(jb,jgrd)323 zflag = idx_bdy( ib_bdy)%flagv(jb,jgrd)323 DO i_bdy = 1, idx_bdy(jbdy)%nblenrim(jgrd) 324 ji = idx_bdy(jbdy)%nbi(i_bdy,jgrd) 325 jj = idx_bdy(jbdy)%nbj(i_bdy,jgrd) 326 zflag = idx_bdy(jbdy)%flagv(i_bdy,jgrd) 324 327 ! 325 328 IF ( ABS( zflag ) == 1. ) THEN ! northern and southern boundaries … … 328 331 zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj-1) ) ) ! 0 if no ice 329 332 ! 330 ! u_ice = u_ice of the adjacent grid point except if this grid point is ice-free (then u_ice = u_oce)333 ! v_ice = v_ice of the adjacent grid point except if this grid point is ice-free (then do not change v_ice) 331 334 v_ice (ji,jj) = v_ice(ji,jj+1) * 0.5_wp * ABS( zflag + 1._wp ) * zmsk1 + & 332 335 & v_ice(ji,jj-1) * 0.5_wp * ABS( zflag - 1._wp ) * zmsk2 + & 333 & v_ oce(ji,jj ) * ( 1._wp - MIN( 1._wp, zmsk1 + zmsk2 ) )336 & v_ice(ji,jj ) * ( 1._wp - MIN( 1._wp, zmsk1 + zmsk2 ) ) 334 337 ELSE ! everywhere else 335 !v_ice(ji,jj) = v_oce(ji,jj)336 338 v_ice(ji,jj) = 0._wp 337 339 ENDIF 338 ! mask ice velocities339 rswitch = MAX( 0.0_wp , SIGN ( 1.0_wp , at_i(ji,jj) - 0.01 ) ) ! 0 if no ice340 v_ice(ji,jj) = rswitch * v_ice(ji,jj)341 340 ! 342 341 END DO 343 CALL lbc_bdy_lnk( v_ice(:,:), 'V', -1., ib_bdy )342 CALL lbc_bdy_lnk( v_ice(:,:), 'V', -1., jbdy ) 344 343 ! 345 344 END SELECT … … 350 349 ! 351 350 END DO 351 ! 352 IF( ln_timing ) CALL timing_stop('bdy_ice_dyn') 352 353 ! 353 354 END SUBROUTINE bdy_ice_dyn … … 359 360 CONTAINS 360 361 SUBROUTINE bdy_ice( kt ) ! Empty routine 362 IMPLICIT NONE 363 INTEGER, INTENT( in ) :: kt 361 364 WRITE(*,*) 'bdy_ice: You should not have seen this print! error?', kt 362 365 END SUBROUTINE bdy_ice
Note: See TracChangeset
for help on using the changeset viewer.