Changeset 11536 for NEMO/trunk/src/OCE/BDY/bdyice.F90
- Timestamp:
- 2019-09-11T15:54:18+02:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/OCE/BDY/bdyice.F90
r11041 r11536 55 55 INTEGER, INTENT(in) :: kt ! Main time step counter 56 56 ! 57 INTEGER :: jbdy ! BDY set index 57 INTEGER :: jbdy, ir ! BDY set index, rim index 58 INTEGER :: ibeg, iend ! length of rim to be treated (rim 0 or rim 1) 59 LOGICAL :: llrim0 ! indicate if rim 0 is treated 60 LOGICAL, DIMENSION(4) :: llsend1, llrecv1 ! indicate how communications are to be carried out 58 61 !!---------------------------------------------------------------------- 59 62 ! controls 60 63 IF( ln_timing ) CALL timing_start('bdy_ice_thd') ! timing 61 64 IF( ln_icediachk ) CALL ice_cons_hsm(0,'bdy_ice_thd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 65 IF( ln_icediachk ) CALL ice_cons2D (0,'bdy_ice_thd', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft) ! conservation 62 66 ! 63 67 CALL ice_var_glo2eqv 64 68 ! 65 DO jbdy = 1, nb_bdy 69 llsend1(:) = .false. ; llrecv1(:) = .false. 70 DO ir = 1, 0, -1 ! treat rim 1 before rim 0 71 IF( ir == 0 ) THEN ; llrim0 = .TRUE. 72 ELSE ; llrim0 = .FALSE. 73 END IF 74 DO jbdy = 1, nb_bdy 75 ! 76 SELECT CASE( cn_ice(jbdy) ) 77 CASE('none') ; CYCLE 78 CASE('frs' ) ; CALL bdy_ice_frs( idx_bdy(jbdy), dta_bdy(jbdy), kt, jbdy, llrim0 ) 79 CASE DEFAULT 80 CALL ctl_stop( 'bdy_ice : unrecognised option for open boundaries for ice fields' ) 81 END SELECT 82 ! 83 END DO 66 84 ! 67 SELECT CASE( cn_ice(jbdy) ) 68 CASE('none') ; CYCLE 69 CASE('frs' ) ; CALL bdy_ice_frs( idx_bdy(jbdy), dta_bdy(jbdy), kt, jbdy ) 70 CASE DEFAULT 71 CALL ctl_stop( 'bdy_ice : unrecognised option for open boundaries for ice fields' ) 72 END SELECT 73 ! 74 END DO 85 ! Update bdy points 86 IF( nn_hls > 1 .AND. ir == 1 ) CYCLE ! at least 2 halos will be corrected -> no need to correct rim 1 before rim 0 87 IF( nn_hls == 1 ) THEN ; llsend1(:) = .false. ; llrecv1(:) = .false. ; END IF 88 DO jbdy = 1, nb_bdy 89 IF( cn_ice(jbdy) == 'frs' ) THEN 90 llsend1(:) = llsend1(:) .OR. lsend_bdyint(jbdy,1,:,ir) ! possibly every direction, T points 91 llrecv1(:) = llrecv1(:) .OR. lrecv_bdyint(jbdy,1,:,ir) ! possibly every direction, T points 92 END IF 93 END DO ! jbdy 94 IF( ANY(llsend1) .OR. ANY(llrecv1) ) THEN ! if need to send/recv in at least one direction 95 ! exchange 3d arrays 96 CALL lbc_lnk_multi( 'bdyice', a_i , 'T', 1., h_i , 'T', 1., h_s , 'T', 1., oa_i, 'T', 1. & 97 & , a_ip, 'T', 1., v_ip, 'T', 1., s_i , 'T', 1., t_su, 'T', 1. & 98 & , v_i , 'T', 1., v_s , 'T', 1., sv_i, 'T', 1. & 99 & , kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 ) 100 ! exchange 4d arrays : third dimension = 1 and then third dimension = jpk 101 CALL lbc_lnk_multi( 'bdyice', t_s , 'T', 1., e_s , 'T', 1., kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 ) 102 CALL lbc_lnk_multi( 'bdyice', t_i , 'T', 1., e_i , 'T', 1., kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 ) 103 END IF 104 END DO ! ir 75 105 ! 76 106 CALL ice_cor( kt , 0 ) ! -- In case categories are out of bounds, do a remapping … … 80 110 ! 81 111 ! controls 112 IF( ln_icectl ) CALL ice_prt ( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' ) ! prints 82 113 IF( ln_icediachk ) CALL ice_cons_hsm(1,'bdy_ice_thd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 83 IF( ln_ice ctl ) CALL ice_prt ( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' ) ! prints114 IF( ln_icediachk ) CALL ice_cons2D (1,'bdy_ice_thd', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft) ! conservation 84 115 IF( ln_timing ) CALL timing_stop ('bdy_ice_thd') ! timing 85 116 ! … … 87 118 88 119 89 SUBROUTINE bdy_ice_frs( idx, dta, kt, jbdy )120 SUBROUTINE bdy_ice_frs( idx, dta, kt, jbdy, llrim0 ) 90 121 !!------------------------------------------------------------------------------ 91 122 !! *** SUBROUTINE bdy_ice_frs *** … … 96 127 !! dimensional baroclinic ocean model with realistic topography. Tellus, 365-382. 97 128 !!------------------------------------------------------------------------------ 98 TYPE(OBC_INDEX), INTENT(in) :: idx ! OBC indices 99 TYPE(OBC_DATA), INTENT(in) :: dta ! OBC external data 100 INTEGER, INTENT(in) :: kt ! main time-step counter 101 INTEGER, INTENT(in) :: jbdy ! BDY set index 129 TYPE(OBC_INDEX), INTENT(in) :: idx ! OBC indices 130 TYPE(OBC_DATA), INTENT(in) :: dta ! OBC external data 131 INTEGER, INTENT(in) :: kt ! main time-step counter 132 INTEGER, INTENT(in) :: jbdy ! BDY set index 133 LOGICAL, INTENT(in) :: llrim0 ! indicate if rim 0 is treated 102 134 ! 103 135 INTEGER :: jpbound ! 0 = incoming ice 104 136 ! ! 1 = outgoing ice 137 INTEGER :: ibeg, iend ! length of rim to be treated (rim 0 or rim 1) 105 138 INTEGER :: i_bdy, jgrd ! dummy loop indices 106 139 INTEGER :: ji, jj, jk, jl, ib, jb 107 140 REAL(wp) :: zwgt, zwgt1 ! local scalar 108 141 REAL(wp) :: ztmelts, zdh 142 REAL(wp), POINTER :: flagu, flagv ! short cuts 109 143 !!------------------------------------------------------------------------------ 110 144 ! 111 145 jgrd = 1 ! Everything is at T-points here 146 IF( llrim0 ) THEN ; ibeg = 1 ; iend = idx%nblenrim0(jgrd) 147 ELSE ; ibeg = idx%nblenrim0(jgrd)+1 ; iend = idx%nblenrim(jgrd) 148 END IF 112 149 ! 113 150 DO jl = 1, jpl 114 DO i_bdy = 1, idx%nblenrim(jgrd)151 DO i_bdy = ibeg, iend 115 152 ji = idx%nbi(i_bdy,jgrd) 116 153 jj = idx%nbj(i_bdy,jgrd) 117 154 zwgt = idx%nbw(i_bdy,jgrd) 118 155 zwgt1 = 1.e0 - idx%nbw(i_bdy,jgrd) 119 a_i(ji,jj,jl) = ( a_i(ji,jj,jl) * zwgt1 + dta%a_i(i_bdy,jl) * zwgt ) * tmask(ji,jj,1) ! Leads fraction 120 h_i(ji,jj,jl) = ( h_i(ji,jj,jl) * zwgt1 + dta%h_i(i_bdy,jl) * zwgt ) * tmask(ji,jj,1) ! Ice depth 121 h_s(ji,jj,jl) = ( h_s(ji,jj,jl) * zwgt1 + dta%h_s(i_bdy,jl) * zwgt ) * tmask(ji,jj,1) ! Snow depth 122 156 a_i (ji,jj, jl) = ( a_i (ji,jj, jl) * zwgt1 + dta%a_i(i_bdy,jl) * zwgt ) * tmask(ji,jj,1) ! Ice concentration 157 h_i (ji,jj, jl) = ( h_i (ji,jj, jl) * zwgt1 + dta%h_i(i_bdy,jl) * zwgt ) * tmask(ji,jj,1) ! Ice depth 158 h_s (ji,jj, jl) = ( h_s (ji,jj, jl) * zwgt1 + dta%h_s(i_bdy,jl) * zwgt ) * tmask(ji,jj,1) ! Snow depth 159 t_i (ji,jj,:,jl) = ( t_i (ji,jj,:,jl) * zwgt1 + dta%t_i(i_bdy,jl) * zwgt ) * tmask(ji,jj,1) ! Ice temperature 160 t_s (ji,jj,:,jl) = ( t_s (ji,jj,:,jl) * zwgt1 + dta%t_s(i_bdy,jl) * zwgt ) * tmask(ji,jj,1) ! Snow temperature 161 t_su(ji,jj, jl) = ( t_su(ji,jj, jl) * zwgt1 + dta%tsu(i_bdy,jl) * zwgt ) * tmask(ji,jj,1) ! Surf temperature 162 s_i (ji,jj, jl) = ( s_i (ji,jj, jl) * zwgt1 + dta%s_i(i_bdy,jl) * zwgt ) * tmask(ji,jj,1) ! Ice salinity 163 a_ip(ji,jj, jl) = ( a_ip(ji,jj, jl) * zwgt1 + dta%aip(i_bdy,jl) * zwgt ) * tmask(ji,jj,1) ! Ice pond concentration 164 h_ip(ji,jj, jl) = ( h_ip(ji,jj, jl) * zwgt1 + dta%hip(i_bdy,jl) * zwgt ) * tmask(ji,jj,1) ! Ice pond depth 165 ! 166 sz_i(ji,jj,:,jl) = s_i(ji,jj,jl) 167 ! 168 ! make sure ponds = 0 if no ponds scheme 169 IF( .NOT.ln_pnd ) THEN 170 a_ip(ji,jj,jl) = 0._wp 171 h_ip(ji,jj,jl) = 0._wp 172 ENDIF 173 ! 123 174 ! ----------------- 124 175 ! Pathological case … … 135 186 h_i(ji,jj,jl) = MIN( hi_max(jl), h_i(ji,jj,jl) + zdh ) 136 187 h_s(ji,jj,jl) = MAX( 0._wp, h_s(ji,jj,jl) - zdh * rhoi / rhos ) 137 188 ! 138 189 ENDDO 139 190 ENDDO 140 CALL lbc_bdy_lnk( 'bdyice', a_i(:,:,:), 'T', 1., jbdy )141 CALL lbc_bdy_lnk( 'bdyice', h_i(:,:,:), 'T', 1., jbdy )142 CALL lbc_bdy_lnk( 'bdyice', h_s(:,:,:), 'T', 1., jbdy )143 191 144 192 DO jl = 1, jpl 145 DO i_bdy = 1, idx%nblenrim(jgrd)193 DO i_bdy = ibeg, iend 146 194 ji = idx%nbi(i_bdy,jgrd) 147 195 jj = idx%nbj(i_bdy,jgrd) 148 196 flagu => idx%flagu(i_bdy,jgrd) 197 flagv => idx%flagv(i_bdy,jgrd) 149 198 ! condition on ice thickness depends on the ice velocity 150 199 ! if velocity is outward (strictly), then ice thickness, volume... must be equal to adjacent values 151 200 jpbound = 0 ; ib = ji ; jb = jj 152 201 ! 153 IF( u_ice(ji ,jj ) < 0. .AND. umask(ji-1,jj ,1) == 0. ) jpbound = 1 ; ib = ji+1 154 IF( u_ice(ji-1,jj ) > 0. .AND. umask(ji ,jj ,1) == 0. ) jpbound = 1 ; ib = ji-1 155 IF( v_ice(ji ,jj ) < 0. .AND. vmask(ji ,jj-1,1) == 0. ) jpbound = 1 ; jb = jj+1 156 IF( v_ice(ji ,jj-1) > 0. .AND. vmask(ji ,jj ,1) == 0. ) jpbound = 1 ; jb = jj-1 202 IF( flagu == 1. ) THEN 203 IF( ji+1 > jpi ) CYCLE 204 IF( u_ice(ji ,jj ) < 0. ) jpbound = 1 ; ib = ji+1 205 END IF 206 IF( flagu == -1. ) THEN 207 IF( ji-1 < 1 ) CYCLE 208 IF( u_ice(ji-1,jj ) < 0. ) jpbound = 1 ; ib = ji-1 209 END IF 210 IF( flagv == 1. ) THEN 211 IF( jj+1 > jpj ) CYCLE 212 IF( v_ice(ji ,jj ) < 0. ) jpbound = 1 ; jb = jj+1 213 END IF 214 IF( flagv == -1. ) THEN 215 IF( jj-1 < 1 ) CYCLE 216 IF( v_ice(ji ,jj-1) < 0. ) jpbound = 1 ; jb = jj-1 217 END IF 157 218 ! 158 219 IF( nn_ice_dta(jbdy) == 0 ) jpbound = 0 ; ib = ji ; jb = jj ! case ice boundaries = initial conditions … … 161 222 IF( a_i(ib,jb,jl) > 0._wp ) THEN ! there is ice at the boundary 162 223 ! 163 a_i(ji,jj,jl) = a_i(ib,jb,jl) ! concentration 164 h_i(ji,jj,jl) = h_i(ib,jb,jl) ! thickness ice 165 h_s(ji,jj,jl) = h_s(ib,jb,jl) ! thickness snw 166 ! 167 SELECT CASE( jpbound ) 168 ! 169 CASE( 0 ) ! velocity is inward 170 ! 171 oa_i(ji,jj, jl) = rn_ice_age(jbdy) * a_i(ji,jj,jl) ! age 172 a_ip(ji,jj, jl) = 0._wp ! pond concentration 173 v_ip(ji,jj, jl) = 0._wp ! pond volume 174 t_su(ji,jj, jl) = rn_ice_tem(jbdy) ! temperature surface 175 t_s (ji,jj,:,jl) = rn_ice_tem(jbdy) ! temperature snw 176 t_i (ji,jj,:,jl) = rn_ice_tem(jbdy) ! temperature ice 177 s_i (ji,jj, jl) = rn_ice_sal(jbdy) ! salinity 178 sz_i(ji,jj,:,jl) = rn_ice_sal(jbdy) ! salinity profile 179 ! 180 CASE( 1 ) ! velocity is outward 181 ! 182 oa_i(ji,jj, jl) = oa_i(ib,jb, jl) ! age 183 a_ip(ji,jj, jl) = a_ip(ib,jb, jl) ! pond concentration 184 v_ip(ji,jj, jl) = v_ip(ib,jb, jl) ! pond volume 185 t_su(ji,jj, jl) = t_su(ib,jb, jl) ! temperature surface 186 t_s (ji,jj,:,jl) = t_s (ib,jb,:,jl) ! temperature snw 187 t_i (ji,jj,:,jl) = t_i (ib,jb,:,jl) ! temperature ice 188 s_i (ji,jj, jl) = s_i (ib,jb, jl) ! salinity 189 sz_i(ji,jj,:,jl) = sz_i(ib,jb,:,jl) ! salinity profile 190 ! 191 END SELECT 224 a_i (ji,jj, jl) = a_i (ib,jb, jl) 225 h_i (ji,jj, jl) = h_i (ib,jb, jl) 226 h_s (ji,jj, jl) = h_s (ib,jb, jl) 227 t_i (ji,jj,:,jl) = t_i (ib,jb,:,jl) 228 t_s (ji,jj,:,jl) = t_s (ib,jb,:,jl) 229 t_su(ji,jj, jl) = t_su(ib,jb, jl) 230 s_i (ji,jj, jl) = s_i (ib,jb, jl) 231 a_ip(ji,jj, jl) = a_ip(ib,jb, jl) 232 h_ip(ji,jj, jl) = h_ip(ib,jb, jl) 233 ! 234 sz_i(ji,jj,:,jl) = sz_i(ib,jb,:,jl) 235 ! 236 ! ice age 237 IF ( jpbound == 0 ) THEN ! velocity is inward 238 oa_i(ji,jj,jl) = rice_age(jbdy) * a_i(ji,jj,jl) 239 ELSEIF( jpbound == 1 ) THEN ! velocity is outward 240 oa_i(ji,jj,jl) = oa_i(ib,jb,jl) 241 ENDIF 192 242 ! 193 243 IF( nn_icesal == 1 ) THEN ! if constant salinity … … 214 264 END DO 215 265 ! 266 ! melt ponds 267 IF( a_i(ji,jj,jl) > epsi10 ) THEN 268 a_ip_frac(ji,jj,jl) = a_ip(ji,jj,jl) / a_i (ji,jj,jl) 269 ELSE 270 a_ip_frac(ji,jj,jl) = 0._wp 271 ENDIF 272 v_ip(ji,jj,jl) = h_ip(ji,jj,jl) * a_ip(ji,jj,jl) 273 ! 216 274 ELSE ! no ice at the boundary 217 275 ! … … 225 283 t_s (ji,jj,:,jl) = rt0 226 284 t_i (ji,jj,:,jl) = rt0 285 286 a_ip_frac(ji,jj,jl) = 0._wp 287 h_ip (ji,jj,jl) = 0._wp 288 a_ip (ji,jj,jl) = 0._wp 289 v_ip (ji,jj,jl) = 0._wp 227 290 228 291 IF( nn_icesal == 1 ) THEN ! if constant salinity … … 246 309 ! 247 310 END DO ! jl 248 249 CALL lbc_bdy_lnk( 'bdyice', a_i (:,:,:) , 'T', 1., jbdy )250 CALL lbc_bdy_lnk( 'bdyice', h_i (:,:,:) , 'T', 1., jbdy )251 CALL lbc_bdy_lnk( 'bdyice', h_s (:,:,:) , 'T', 1., jbdy )252 CALL lbc_bdy_lnk( 'bdyice', oa_i(:,:,:) , 'T', 1., jbdy )253 CALL lbc_bdy_lnk( 'bdyice', a_ip(:,:,:) , 'T', 1., jbdy )254 CALL lbc_bdy_lnk( 'bdyice', v_ip(:,:,:) , 'T', 1., jbdy )255 CALL lbc_bdy_lnk( 'bdyice', s_i (:,:,:) , 'T', 1., jbdy )256 CALL lbc_bdy_lnk( 'bdyice', t_su(:,:,:) , 'T', 1., jbdy )257 CALL lbc_bdy_lnk( 'bdyice', v_i (:,:,:) , 'T', 1., jbdy )258 CALL lbc_bdy_lnk( 'bdyice', v_s (:,:,:) , 'T', 1., jbdy )259 CALL lbc_bdy_lnk( 'bdyice', sv_i(:,:,:) , 'T', 1., jbdy )260 CALL lbc_bdy_lnk( 'bdyice', t_s (:,:,:,:), 'T', 1., jbdy )261 CALL lbc_bdy_lnk( 'bdyice', e_s (:,:,:,:), 'T', 1., jbdy )262 CALL lbc_bdy_lnk( 'bdyice', t_i (:,:,:,:), 'T', 1., jbdy )263 CALL lbc_bdy_lnk( 'bdyice', e_i (:,:,:,:), 'T', 1., jbdy )264 311 ! 265 312 END SUBROUTINE bdy_ice_frs … … 279 326 CHARACTER(len=1), INTENT(in) :: cd_type ! nature of velocity grid-points 280 327 ! 281 INTEGER :: i_bdy, jgrd ! dummy loop indices 282 INTEGER :: ji, jj ! local scalar 283 INTEGER :: jbdy ! BDY set index 328 INTEGER :: i_bdy, jgrd ! dummy loop indices 329 INTEGER :: ji, jj ! local scalar 330 INTEGER :: jbdy, ir ! BDY set index, rim index 331 INTEGER :: ibeg, iend ! length of rim to be treated (rim 0 or rim 1) 284 332 REAL(wp) :: zmsk1, zmsk2, zflag 333 LOGICAL, DIMENSION(4) :: llsend2, llrecv2, llsend3, llrecv3 ! indicate how communications are to be carried out 285 334 !!------------------------------------------------------------------------------ 286 335 IF( ln_timing ) CALL timing_start('bdy_ice_dyn') 287 336 ! 288 DO jbdy=1, nb_bdy 337 llsend2(:) = .false. ; llrecv2(:) = .false. 338 llsend3(:) = .false. ; llrecv3(:) = .false. 339 DO ir = 1, 0, -1 340 DO jbdy = 1, nb_bdy 341 ! 342 SELECT CASE( cn_ice(jbdy) ) 343 ! 344 CASE('none') 345 CYCLE 346 ! 347 CASE('frs') 348 ! 349 IF( nn_ice_dta(jbdy) == 0 ) CYCLE ! case ice boundaries = initial conditions 350 ! ! do not change ice velocity (it is only computed by rheology) 351 SELECT CASE ( cd_type ) 352 ! 353 CASE ( 'U' ) 354 jgrd = 2 ! u velocity 355 IF( ir == 0 ) THEN ; ibeg = 1 ; iend = idx_bdy(jbdy)%nblenrim0(jgrd) 356 ELSE ; ibeg = idx_bdy(jbdy)%nblenrim0(jgrd)+1 ; iend = idx_bdy(jbdy)%nblenrim(jgrd) 357 END IF 358 DO i_bdy = ibeg, iend 359 ji = idx_bdy(jbdy)%nbi(i_bdy,jgrd) 360 jj = idx_bdy(jbdy)%nbj(i_bdy,jgrd) 361 zflag = idx_bdy(jbdy)%flagu(i_bdy,jgrd) 362 ! i-1 i i | ! i i i+1 | ! i i i+1 | 363 ! > ice > | ! o > ice | ! o > o | 364 ! => set at u_ice(i-1) ! => set to O ! => unchanged 365 IF( zflag == -1. .AND. ji > 1 .AND. ji < jpi ) THEN 366 IF ( vt_i(ji ,jj) > 0. ) THEN ; u_ice(ji,jj) = u_ice(ji-1,jj) 367 ELSEIF( vt_i(ji+1,jj) > 0. ) THEN ; u_ice(ji,jj) = 0._wp 368 END IF 369 END IF 370 ! | i i+1 i+1 ! | i i i+1 ! | i i i+1 371 ! | > ice > ! | ice > o ! | o > o 372 ! => set at u_ice(i+1) ! => set to O ! => unchanged 373 IF( zflag == 1. .AND. ji+1 < jpi+1 ) THEN 374 IF ( vt_i(ji+1,jj) > 0. ) THEN ; u_ice(ji,jj) = u_ice(ji+1,jj) 375 ELSEIF( vt_i(ji ,jj) > 0. ) THEN ; u_ice(ji,jj) = 0._wp 376 END IF 377 END IF 378 ! 379 IF( zflag == 0. ) u_ice(ji,jj) = 0._wp ! u_ice = 0 if north/south bdy 380 ! 381 END DO 382 ! 383 CASE ( 'V' ) 384 jgrd = 3 ! v velocity 385 IF( ir == 0 ) THEN ; ibeg = 1 ; iend = idx_bdy(jbdy)%nblenrim0(jgrd) 386 ELSE ; ibeg = idx_bdy(jbdy)%nblenrim0(jgrd)+1 ; iend = idx_bdy(jbdy)%nblenrim(jgrd) 387 END IF 388 DO i_bdy = ibeg, iend 389 ji = idx_bdy(jbdy)%nbi(i_bdy,jgrd) 390 jj = idx_bdy(jbdy)%nbj(i_bdy,jgrd) 391 zflag = idx_bdy(jbdy)%flagv(i_bdy,jgrd) 392 ! ! ice (jj+1) ! o (jj+1) 393 ! ^ (jj ) ! ^ (jj ) ! ^ (jj ) 394 ! ice (jj ) ! o (jj ) ! o (jj ) 395 ! ^ (jj-1) ! ! 396 ! => set to u_ice(jj-1) ! => set to 0 ! => unchanged 397 IF( zflag == -1. .AND. jj > 1 .AND. jj < jpj ) THEN 398 IF ( vt_i(ji,jj ) > 0. ) THEN ; v_ice(ji,jj) = v_ice(ji,jj-1) 399 ELSEIF( vt_i(ji,jj+1) > 0. ) THEN ; v_ice(ji,jj) = 0._wp 400 END IF 401 END IF 402 ! ^ (jj+1) ! ! 403 ! ice (jj+1) ! o (jj+1) ! o (jj+1) 404 ! ^ (jj ) ! ^ (jj ) ! ^ (jj ) 405 ! ________________ ! ____ice___(jj )_ ! _____o____(jj ) 406 ! => set to u_ice(jj+1) ! => set to 0 ! => unchanged 407 IF( zflag == 1. .AND. jj < jpj ) THEN 408 IF ( vt_i(ji,jj+1) > 0. ) THEN ; v_ice(ji,jj) = v_ice(ji,jj+1) 409 ELSEIF( vt_i(ji,jj ) > 0. ) THEN ; v_ice(ji,jj) = 0._wp 410 END IF 411 END IF 412 ! 413 IF( zflag == 0. ) v_ice(ji,jj) = 0._wp ! v_ice = 0 if west/east bdy 414 ! 415 END DO 416 ! 417 END SELECT 418 ! 419 CASE DEFAULT 420 CALL ctl_stop( 'bdy_ice_dyn : unrecognised option for open boundaries for ice fields' ) 421 END SELECT 422 ! 423 END DO ! jbdy 289 424 ! 290 SELECT CASE( cn_ice(jbdy) ) 291 ! 292 CASE('none') 293 CYCLE 294 ! 295 CASE('frs') 296 ! 297 IF( nn_ice_dta(jbdy) == 0 ) CYCLE ! case ice boundaries = initial conditions 298 ! ! do not change ice velocity (it is only computed by rheology) 299 SELECT CASE ( cd_type ) 300 ! 301 CASE ( 'U' ) 302 jgrd = 2 ! u velocity 303 DO i_bdy = 1, idx_bdy(jbdy)%nblenrim(jgrd) 304 ji = idx_bdy(jbdy)%nbi(i_bdy,jgrd) 305 jj = idx_bdy(jbdy)%nbj(i_bdy,jgrd) 306 zflag = idx_bdy(jbdy)%flagu(i_bdy,jgrd) 307 ! 308 IF ( ABS( zflag ) == 1. ) THEN ! eastern and western boundaries 309 ! one of the two zmsk is always 0 (because of zflag) 310 zmsk1 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji+1,jj) ) ) ! 0 if no ice 311 zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj) ) ) ! 0 if no ice 312 ! 313 ! u_ice = u_ice of the adjacent grid point except if this grid point is ice-free (then do not change u_ice) 314 u_ice (ji,jj) = u_ice(ji+1,jj) * 0.5_wp * ABS( zflag + 1._wp ) * zmsk1 + & 315 & u_ice(ji-1,jj) * 0.5_wp * ABS( zflag - 1._wp ) * zmsk2 + & 316 & u_ice(ji ,jj) * ( 1._wp - MIN( 1._wp, zmsk1 + zmsk2 ) ) 317 ELSE ! everywhere else 318 u_ice(ji,jj) = 0._wp 319 ENDIF 320 ! 321 END DO 322 CALL lbc_bdy_lnk( 'bdyice', u_ice(:,:), 'U', -1., jbdy ) 323 ! 324 CASE ( 'V' ) 325 jgrd = 3 ! v velocity 326 DO i_bdy = 1, idx_bdy(jbdy)%nblenrim(jgrd) 327 ji = idx_bdy(jbdy)%nbi(i_bdy,jgrd) 328 jj = idx_bdy(jbdy)%nbj(i_bdy,jgrd) 329 zflag = idx_bdy(jbdy)%flagv(i_bdy,jgrd) 330 ! 331 IF ( ABS( zflag ) == 1. ) THEN ! northern and southern boundaries 332 ! one of the two zmsk is always 0 (because of zflag) 333 zmsk1 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj+1) ) ) ! 0 if no ice 334 zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj) ) ) ! 0 if no ice 335 ! 336 ! v_ice = v_ice of the adjacent grid point except if this grid point is ice-free (then do not change v_ice) 337 v_ice (ji,jj) = v_ice(ji,jj+1) * 0.5_wp * ABS( zflag + 1._wp ) * zmsk1 + & 338 & v_ice(ji,jj-1) * 0.5_wp * ABS( zflag - 1._wp ) * zmsk2 + & 339 & v_ice(ji,jj ) * ( 1._wp - MIN( 1._wp, zmsk1 + zmsk2 ) ) 340 ELSE ! everywhere else 341 v_ice(ji,jj) = 0._wp 342 ENDIF 343 ! 344 END DO 345 CALL lbc_bdy_lnk( 'bdyice', v_ice(:,:), 'V', -1., jbdy ) 346 ! 347 END SELECT 348 ! 349 CASE DEFAULT 350 CALL ctl_stop( 'bdy_ice_dyn : unrecognised option for open boundaries for ice fields' ) 425 SELECT CASE ( cd_type ) 426 CASE ( 'U' ) 427 IF( nn_hls > 1 .AND. ir == 1 ) CYCLE ! at least 2 halos will be corrected -> no need to correct rim 1 before rim 0 428 IF( nn_hls == 1 ) THEN ; llsend2(:) = .false. ; llrecv2(:) = .false. ; END IF 429 DO jbdy = 1, nb_bdy 430 IF( cn_ice(jbdy) == 'frs' .AND. nn_ice_dta(jbdy) /= 0 ) THEN 431 llsend2(:) = llsend2(:) .OR. lsend_bdyint(jbdy,2,:,ir) ! possibly every direction, U points 432 llsend2(1) = llsend2(1) .OR. lsend_bdyext(jbdy,2,1,ir) ! neighbour might search point towards its west bdy 433 llrecv2(:) = llrecv2(:) .OR. lrecv_bdyint(jbdy,2,:,ir) ! possibly every direction, U points 434 llrecv2(2) = llrecv2(2) .OR. lrecv_bdyext(jbdy,2,2,ir) ! might search point towards east bdy 435 END IF 436 END DO 437 IF( ANY(llsend2) .OR. ANY(llrecv2) ) THEN ! if need to send/recv in at least one direction 438 CALL lbc_lnk( 'bdyice', u_ice, 'U', -1., kfillmode=jpfillnothing ,lsend=llsend2, lrecv=llrecv2 ) 439 END IF 440 CASE ( 'V' ) 441 IF( nn_hls > 1 .AND. ir == 1 ) CYCLE ! at least 2 halos will be corrected -> no need to correct rim 1 before rim 0 442 IF( nn_hls == 1 ) THEN ; llsend3(:) = .false. ; llrecv3(:) = .false. ; END IF 443 DO jbdy = 1, nb_bdy 444 IF( cn_ice(jbdy) == 'frs' .AND. nn_ice_dta(jbdy) /= 0 ) THEN 445 llsend3(:) = llsend3(:) .OR. lsend_bdyint(jbdy,3,:,ir) ! possibly every direction, V points 446 llsend3(3) = llsend3(3) .OR. lsend_bdyext(jbdy,3,3,ir) ! neighbour might search point towards its south bdy 447 llrecv3(:) = llrecv3(:) .OR. lrecv_bdyint(jbdy,3,:,ir) ! possibly every direction, V points 448 llrecv3(4) = llrecv3(4) .OR. lrecv_bdyext(jbdy,3,4,ir) ! might search point towards north bdy 449 END IF 450 END DO 451 IF( ANY(llsend3) .OR. ANY(llrecv3) ) THEN ! if need to send/recv in at least one direction 452 CALL lbc_lnk( 'bdyice', v_ice, 'V', -1., kfillmode=jpfillnothing ,lsend=llsend3, lrecv=llrecv3 ) 453 END IF 351 454 END SELECT 352 ! 353 END DO 455 END DO ! ir 354 456 ! 355 457 IF( ln_timing ) CALL timing_stop('bdy_ice_dyn')
Note: See TracChangeset
for help on using the changeset viewer.