Changeset 6416 for trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
- Timestamp:
- 2016-04-01T14:22:17+02:00 (8 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90 (modified) (11 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
r6403 r6416 75 75 INTEGER :: ii, ij, iter ! - - 76 76 REAL(wp) :: ztmelts, zdv, zfrazb, zweight, zde ! local scalars 77 REAL(wp) :: zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf , zhicol_new! - -77 REAL(wp) :: zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf ! - - 78 78 REAL(wp) :: ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit ! - - 79 LOGICAL :: iterate_frazil ! iterate frazil ice collection thickness80 79 CHARACTER (len = 15) :: fieldid 81 80 … … 108 107 REAL(wp), POINTER, DIMENSION(:,:) :: zsmv_i_1d ! 1-D version of smv_i 109 108 110 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze_i_1d 111 112 REAL(wp), POINTER, DIMENSION(:,:) :: zvrel 113 114 REAL(wp) :: zcai = 1.4e-3_wp 109 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze_i_1d !: 1-D version of e_i 110 111 REAL(wp), POINTER, DIMENSION(:,:) :: zvrel ! relative ice / frazil velocity 112 113 REAL(wp) :: zcai = 1.4e-3_wp ! ice-air drag (clem: should be dependent on coupling/forcing used) 115 114 !!-----------------------------------------------------------------------! 116 115 … … 143 142 !------------------------------------------------------------------------------! 144 143 ! hicol is the thickness of new ice formed in open water 145 ! hicol can be either prescribed (frazswi = 0) 146 ! or computed (frazswi = 1) 144 ! hicol can be either prescribed (frazswi = 0) or computed (frazswi = 1) 147 145 ! Frazil ice forms in open water, is transported by wind 148 146 ! accumulates at the edge of the consolidated ice edge … … 155 153 zvrel(:,:) = 0._wp 156 154 157 ! Default new ice thickness 158 hicol(:,:) = rn_hnewice 155 ! Default new ice thickness 156 WHERE( qlead(:,:) < 0._wp ) ; hicol = rn_hnewice 157 ELSEWHERE ; hicol = 0._wp 158 END WHERE 159 159 160 160 IF( ln_frazil ) THEN … … 182 182 & + vtau_ice(ji ,jj ) * vmask(ji ,jj ,1) ) * 0.5_wp 183 183 ! Square root of wind stress 184 ztenagm = SQRT( SQRT( ztaux **2 + ztauy**2) )184 ztenagm = SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) ) 185 185 186 186 !--------------------- … … 205 205 zvrel2 = MAX( ( zvfrx - zvgx ) * ( zvfrx - zvgx ) & 206 206 & + ( zvfry - zvgy ) * ( zvfry - zvgy ) , 0.15 * 0.15 ) 207 zvrel(ji,jj) 207 zvrel(ji,jj) = SQRT( zvrel2 ) 208 208 209 209 !--------------------- 210 210 ! Iterative procedure 211 211 !--------------------- 212 hicol(ji,jj) = zhicrit + 0.1 213 hicol(ji,jj) = zhicrit + hicol(ji,jj) & 214 & / ( hicol(ji,jj) * hicol(ji,jj) - zhicrit * zhicrit ) * ztwogp * zvrel2 215 216 !!gm better coding: above: hicol(ji,jj) * hicol(ji,jj) = (zhicrit + 0.1)*(zhicrit + 0.1) 217 !!gm = zhicrit**2 + 0.2*zhicrit +0.01 218 !!gm therefore the 2 lines with hicol can be replaced by 1 line: 219 !!gm hicol(ji,jj) = zhicrit + (zhicrit + 0.1) / ( 0.2 * zhicrit + 0.01 ) * ztwogp * zvrel2 220 !!gm further more (zhicrit + 0.1)/(0.2 * zhicrit + 0.01 )*ztwogp can be computed one for all outside the DO loop 212 hicol(ji,jj) = zhicrit + ( zhicrit + 0.1 ) & 213 & / ( ( zhicrit + 0.1 ) * ( zhicrit + 0.1 ) - zhicrit * zhicrit ) * ztwogp * zvrel2 221 214 222 215 iter = 1 223 iterate_frazil = .true. 224 225 DO WHILE ( iter < 100 .AND. iterate_frazil ) 226 zf = ( hicol(ji,jj) - zhicrit ) * ( hicol(ji,jj)**2 - zhicrit**2 ) & 227 - hicol(ji,jj) * zhicrit * ztwogp * zvrel2 228 zfp = ( hicol(ji,jj) - zhicrit ) * ( 3.0*hicol(ji,jj) + zhicrit ) & 229 - zhicrit * ztwogp * zvrel2 230 zhicol_new = hicol(ji,jj) - zf/zfp 231 hicol(ji,jj) = zhicol_new 232 216 DO WHILE ( iter < 20 ) 217 zf = ( hicol(ji,jj) - zhicrit ) * ( hicol(ji,jj) * hicol(ji,jj) - zhicrit * zhicrit ) - & 218 & hicol(ji,jj) * zhicrit * ztwogp * zvrel2 219 zfp = ( hicol(ji,jj) - zhicrit ) * ( 3.0 * hicol(ji,jj) + zhicrit ) - zhicrit * ztwogp * zvrel2 220 221 hicol(ji,jj) = hicol(ji,jj) - zf/zfp 233 222 iter = iter + 1 234 235 END DO ! do while 223 END DO 236 224 237 225 ENDIF ! end of selection of pixels where ice forms 238 226 239 END DO ! loop on ji ends240 END DO ! loop on jj ends241 !242 CALL lbc_lnk( zvrel(:,:), 'T', 1. )243 CALL lbc_lnk( hicol(:,:), 'T', 1. )227 END DO 228 END DO 229 ! 230 CALL lbc_lnk( zvrel(:,:), 'T', 1. ) 231 CALL lbc_lnk( hicol(:,:), 'T', 1. ) 244 232 245 233 ENDIF ! End of computation of frazil ice collection thickness … … 282 270 ! Move from 2-D to 1-D vectors 283 271 !------------------------------ 284 ! If ocean gains heat do nothing 285 ! 0therwise compute new ice formation 272 ! If ocean gains heat do nothing. Otherwise compute new ice formation 286 273 287 274 IF ( nbpac > 0 ) THEN … … 317 304 zv_b(1:nbpac,:) = zv_i_1d(1:nbpac,:) 318 305 za_b(1:nbpac,:) = za_i_1d(1:nbpac,:) 306 319 307 !---------------------- 320 308 ! Thickness of new ice 321 309 !---------------------- 322 DO ji = 1, nbpac 323 zh_newice(ji) = rn_hnewice 324 END DO 325 IF( ln_frazil ) zh_newice(1:nbpac) = hicol_1d(1:nbpac) 310 zh_newice(1:nbpac) = hicol_1d(1:nbpac) 326 311 327 312 !---------------------- … … 347 332 DO ji = 1, nbpac 348 333 ztmelts = - tmut * zs_newice(ji) + rt0 ! Melting point (K) 349 ze_newice(ji) = rhoic * ( cpic * ( ztmelts - t_bo_1d(ji) ) &334 ze_newice(ji) = rhoic * ( cpic * ( ztmelts - t_bo_1d(ji) ) & 350 335 & + lfus * ( 1.0 - ( ztmelts - rt0 ) / MIN( t_bo_1d(ji) - rt0, -epsi10 ) ) & 351 336 & - rcp * ( ztmelts - rt0 ) ) … … 385 370 ! salt flux 386 371 sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice(ji) * rhoic * zs_newice(ji) * r1_rdtice 387 372 END DO 373 374 zv_frazb(:) = 0._wp 375 IF( ln_frazil ) THEN 388 376 ! A fraction zfrazb of frazil ice is accreted at the ice bottom 389 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp , - zat_i_1d(ji) ) ) 390 zfrazb = rswitch * ( TANH ( rn_Cfrazb * ( zvrel_1d(ji) - rn_vfrazb ) ) + 1.0 ) * 0.5 * rn_maxfrazb 391 zv_frazb(ji) = zfrazb * zv_newice(ji) 392 zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji) 393 END DO 394 377 DO ji = 1, nbpac 378 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp , - zat_i_1d(ji) ) ) 379 zfrazb = rswitch * ( TANH( rn_Cfrazb * ( zvrel_1d(ji) - rn_vfrazb ) ) + 1.0 ) * 0.5 * rn_maxfrazb 380 zv_frazb(ji) = zfrazb * zv_newice(ji) 381 zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji) 382 END DO 383 END IF 384 395 385 !----------------- 396 386 ! Area of new ice … … 444 434 jl = jcat(ji) 445 435 rswitch = MAX( 0._wp, SIGN( 1._wp , zv_i_1d(ji,jl) - epsi20 ) ) 446 ze_i_1d(ji,jk,jl) = zswinew(ji) * ze_newice(ji) + 436 ze_i_1d(ji,jk,jl) = zswinew(ji) * ze_newice(ji) + & 447 437 & ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_1d(ji,jk,jl) * zv_b(ji,jl) ) & 448 438 & * rswitch / MAX( zv_i_1d(ji,jl), epsi20 )
Note: See TracChangeset
for help on using the changeset viewer.
