Changeset 15334 for NEMO/trunk/src/ICE/icethd.F90
 Timestamp:
 20211005T23:18:34+02:00 (14 months ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

NEMO/trunk/src/ICE/icethd.F90
r14997 r15334 20 20 USE sbc_oce , ONLY : sss_m, sst_m, e3t_m, utau, vtau, ssu_m, ssv_m, frq_m, sprecip, ln_cpl 21 21 USE sbc_ice , ONLY : qsr_oce, qns_oce, qemp_oce, qsr_ice, qns_ice, dqns_ice, evap_ice, qprec_ice, qevap_ice, & 22 & qml_ice, qcn_ice, qtr_ice_top 22 & qml_ice, qcn_ice, qtr_ice_top, utau_ice, vtau_ice 23 23 USE ice1D ! seaice: thermodynamics variables 24 24 USE icethd_zdf ! seaice: vertical heat diffusion … … 97 97 REAL(wp), DIMENSION(jpi,jpj) :: zu_io, zv_io, zfric, zvel ! iceocean velocity (m/s) and frictional velocity (m2/s2) 98 98 ! 99 ! for collection thickness 100 INTEGER :: iter !   101 REAL(wp) :: zvfrx, zvgx, ztaux, zf !   102 REAL(wp) :: ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, ztwogp !   103 REAL(wp), PARAMETER :: zcai = 1.4e3_wp ! iceair drag (clem: should be dependent on coupling/forcing used) 104 REAL(wp), PARAMETER :: zhicrit = 0.04_wp ! frazil ice thickness 105 REAL(wp), PARAMETER :: zsqcd = 1.0_wp / SQRT( 1.3_wp * zcai ) ! 1/SQRT(airdensity*drag) 106 REAL(wp), PARAMETER :: zgamafr = 0.03_wp 99 107 !! 100 108 ! controls … … 128 136 & ( v_ice(ji,jj1) + v_ice(ji,jj) ) * ( v_ice(ji,jj1) + v_ice(ji,jj) ) ) 129 137 END_2D 130 ELSE ! if no ice dynamics => trans mitdirectly the atmospheric stress to the ocean138 ELSE ! if no ice dynamics => transfer directly the atmospheric stress to the ocean 131 139 DO_2D( 0, 0, 0, 0 ) 132 140 zfric(ji,jj) = r1_rho0 * SQRT( 0.5_wp * & … … 219 227 ENDIF 220 228 229 230 !! 231 ! Collection thickness of ice formed in leads and polynyas 232 !! 233 ! ht_i_new is the thickness of new ice formed in open water 234 ! ht_i_new can be either prescribed (ln_frazil=F) or computed (ln_frazil=T) 235 ! Frazil ice forms in open water, is transported by wind, accumulates at the edge of the consolidated ice edge 236 ! where it forms aggregates of a specific thickness called collection thickness. 237 ! 238 fraz_frac(:,:) = 0._wp 239 ! 240 ! Default new ice thickness 241 WHERE( qlead(:,:) < 0._wp ) ! cooling 242 ht_i_new(:,:) = rn_hinew 243 ELSEWHERE 244 ht_i_new(:,:) = 0._wp 245 END WHERE 246 247 IF( ln_frazil ) THEN 248 ztwogp = 2._wp * rho0 / ( grav * 0.3_wp * ( rho0  rhoi ) ) ! reduced grav 249 ! 250 DO_2D( 0, 0, 0, 0 ) 251 IF ( qlead(ji,jj) < 0._wp ) THEN ! cooling 252 !  Wind stress  ! 253 ztaux = ( utau_ice(ji1,jj ) * umask(ji1,jj ,1) & 254 & + utau_ice(ji ,jj ) * umask(ji ,jj ,1) ) * 0.5_wp 255 ztauy = ( vtau_ice(ji ,jj1) * vmask(ji ,jj1,1) & 256 & + vtau_ice(ji ,jj ) * vmask(ji ,jj ,1) ) * 0.5_wp 257 ! Square root of wind stress 258 ztenagm = SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) ) 259 260 !  Frazil ice velocity  ! 261 rswitch = MAX( 0._wp, SIGN( 1._wp , ztenagm  epsi10 ) ) 262 zvfrx = rswitch * zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 ) 263 zvfry = rswitch * zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 ) 264 265 !  Pack ice velocity  ! 266 zvgx = ( u_ice(ji1,jj ) * umask(ji1,jj ,1) + u_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp 267 zvgy = ( v_ice(ji ,jj1) * vmask(ji ,jj1,1) + v_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp 268 269 !  Relative frazil/pack ice velocity  ! 270 rswitch = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj)  epsi10 ) ) 271 zvrel2 = MAX( ( zvfrx  zvgx ) * ( zvfrx  zvgx ) & 272 & + ( zvfry  zvgy ) * ( zvfry  zvgy ) , 0.15_wp * 0.15_wp ) * rswitch 273 274 !!clem 275 fraz_frac(ji,jj) = rswitch * ( TANH( rn_Cfraz * ( SQRT(zvrel2)  rn_vfraz ) ) + 1._wp ) * 0.5_wp * rn_maxfraz 276 !!clem 277 278 !  new ice thickness (iterative loop)  ! 279 ht_i_new(ji,jj) = zhicrit + ( zhicrit + 0.1_wp ) & 280 & / ( ( zhicrit + 0.1_wp ) * ( zhicrit + 0.1_wp )  zhicrit * zhicrit ) * ztwogp * zvrel2 281 iter = 1 282 DO WHILE ( iter < 20 ) 283 zf = ( ht_i_new(ji,jj)  zhicrit ) * ( ht_i_new(ji,jj) * ht_i_new(ji,jj)  zhicrit * zhicrit )  & 284 & ht_i_new(ji,jj) * zhicrit * ztwogp * zvrel2 285 zfp = ( ht_i_new(ji,jj)  zhicrit ) * ( 3.0_wp * ht_i_new(ji,jj) + zhicrit )  zhicrit * ztwogp * zvrel2 286 287 ht_i_new(ji,jj) = ht_i_new(ji,jj)  zf / MAX( zfp, epsi20 ) 288 iter = iter + 1 289 END DO 290 ! 291 ! bound ht_i_new (though I don't see why it should be necessary) 292 ht_i_new(ji,jj) = MAX( 0.01_wp, MIN( ht_i_new(ji,jj), hi_max(jpl) ) ) 293 ! 294 ELSE 295 ht_i_new(ji,jj) = 0._wp 296 ENDIF 297 ! 298 END_2D 299 ! 300 CALL lbc_lnk( 'icethd', fraz_frac, 'T', 1.0_wp, ht_i_new, 'T', 1.0_wp ) 301 302 ENDIF 303 221 304 !! 222 305 ! Thermodynamic computation (only on grid points covered by ice) => loop over ice categories … … 268 351 ! 269 352 IF ( ln_pnd .AND. ln_icedH ) & 270 & CALL ice_thd_pnd !  Melt ponds 353 & CALL ice_thd_pnd !  Melt ponds  ! 271 354 ! 272 355 IF( jpl > 1 ) CALL ice_itd_rem( kt ) !  Transport ice between thickness categories  ! … … 276 359 CALL ice_cor( kt , 2 ) !  Corrections  ! 277 360 ! 278 oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rDt_ice ! ice natural aging incrementation 361 oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rDt_ice !  Ice natural aging incrementation 362 ! 363 DO_2D( 0, 0, 0, 0 ) !  Ice velocity corrections 364 IF( at_i(ji,jj) == 0._wp ) THEN ! if ice has melted 365 IF( at_i(ji+1,jj) == 0._wp ) u_ice(ji ,jj) = 0._wp ! right side 366 IF( at_i(ji1,jj) == 0._wp ) u_ice(ji1,jj) = 0._wp ! left side 367 IF( at_i(ji,jj+1) == 0._wp ) v_ice(ji,jj ) = 0._wp ! upper side 368 IF( at_i(ji,jj1) == 0._wp ) v_ice(ji,jj1) = 0._wp ! bottom side 369 ENDIF 370 END_2D 371 CALL lbc_lnk( 'icecor', u_ice, 'U', 1.0_wp, v_ice, 'V', 1.0_wp ) 279 372 ! 280 373 ! convergence tests
Note: See TracChangeset
for help on using the changeset viewer.