Changeset 5208 for branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
- Timestamp:
- 2015-04-13T15:08:59+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
r4688 r5208 37 37 PUBLIC lim_trp ! called by ice_step 38 38 39 REAL(wp) :: epsi10 = 1.e-10_wp40 REAL(wp) :: epsi20 = 1.e-20_wp41 42 39 !! * Substitution 43 40 # include "vectopt_loop_substitute.h90" … … 63 60 INTEGER, INTENT(in) :: kt ! number of iteration 64 61 ! 65 INTEGER :: ji, jj, jk, jl, layer! dummy loop indices62 INTEGER :: ji, jj, jk, jl, jn ! dummy loop indices 66 63 INTEGER :: initad ! number of sub-timestep for the advection 67 INTEGER :: ierr ! error status 68 REAL(wp) :: zindb , zindsn , zindic, zindh, zinda ! local scalar 69 REAL(wp) :: zcfl , zusnit ! - - 70 REAL(wp) :: zsal , zage ! - - 64 REAL(wp) :: zcfl , zusnit ! - - 71 65 ! 72 66 REAL(wp), POINTER, DIMENSION(:,:) :: zui_u, zvi_v, zsm, zs0at, zs0ow 73 67 REAL(wp), POINTER, DIMENSION(:,:,:) :: zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi 74 68 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zs0e 75 ! mass and salt flux (clem) 76 REAL(wp), POINTER, DIMENSION(:,:,:) :: zviold, zvsold ! old ice volume... 77 REAL(wp), POINTER, DIMENSION(:,:,:) :: zaiold, zhimax ! old ice concentration and thickness 78 REAL(wp), POINTER, DIMENSION(:,:) :: zeiold, zesold ! old enthalpies 69 REAL(wp), POINTER, DIMENSION(:,:,:) :: zviold, zvsold ! old ice volume... 70 REAL(wp), POINTER, DIMENSION(:,:,:) :: zaiold, zhimax ! old ice concentration and thickness 71 REAL(wp), POINTER, DIMENSION(:,:) :: zeiold, zesold ! old enthalpies 79 72 REAL(wp) :: zdv, zda, zvi, zvs, zsmv, zes, zei 80 73 ! … … 85 78 CALL wrk_alloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow, zeiold, zesold ) 86 79 CALL wrk_alloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi ) 87 CALL wrk_alloc( jpi, jpj, jkmax, jpl, zs0e )80 CALL wrk_alloc( jpi, jpj, nlay_i+1, jpl, zs0e ) 88 81 89 82 CALL wrk_alloc( jpi, jpj, jpl, zaiold, zhimax, zviold, zvsold ) ! clem … … 167 160 168 161 IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN !== odd ice time step: adv_x then adv_y ==! 169 DO j k= 1,initad162 DO jn = 1,initad 170 163 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0ow (:,:), sxopw(:,:), & !--- ice open water area 171 164 & sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) ) … … 197 190 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), & 198 191 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 199 DO layer= 1, nlay_i !--- ice heat contents ---200 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0e(:,:, layer,jl), sxe (:,:,layer,jl), &201 & sxxe(:,:, layer,jl), sye (:,:,layer,jl), &202 & syye(:,:, layer,jl), sxye(:,:,layer,jl) )203 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0e(:,:, layer,jl), sxe (:,:,layer,jl), &204 & sxxe(:,:, layer,jl), sye (:,:,layer,jl), &205 & syye(:,:, layer,jl), sxye(:,:,layer,jl) )192 DO jk = 1, nlay_i !--- ice heat contents --- 193 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl), & 194 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), & 195 & syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 196 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl), & 197 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), & 198 & syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 206 199 END DO 207 200 END DO 208 201 END DO 209 202 ELSE 210 DO j k= 1, initad203 DO jn = 1, initad 211 204 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0ow (:,:), sxopw(:,:), & !--- ice open water area 212 205 & sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) ) … … 239 232 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), & 240 233 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 241 DO layer= 1, nlay_i !--- ice heat contents ---242 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0e(:,:, layer,jl), sxe (:,:,layer,jl), &243 & sxxe(:,:, layer,jl), sye (:,:,layer,jl), &244 & syye(:,:, layer,jl), sxye(:,:,layer,jl) )245 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0e(:,:, layer,jl), sxe (:,:,layer,jl), &246 & sxxe(:,:, layer,jl), sye (:,:,layer,jl), &247 & syye(:,:, layer,jl), sxye(:,:,layer,jl) )234 DO jk = 1, nlay_i !--- ice heat contents --- 235 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl), & 236 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), & 237 & syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 238 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl), & 239 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), & 240 & syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 248 241 END DO 249 242 END DO … … 341 334 DO jj = 1, jpj 342 335 DO ji = 1, jpi 343 zindb= MAX( 0._wp , SIGN( 1._wp, zs0a(ji,jj,jl) - epsi10 ) )336 rswitch = MAX( 0._wp , SIGN( 1._wp, zs0a(ji,jj,jl) - epsi10 ) ) 344 337 345 338 zvi = zs0ice(ji,jj,jl) … … 349 342 ! 350 343 ! Remove very small areas 351 v_s(ji,jj,jl) = zindb* zs0sn (ji,jj,jl)352 v_i(ji,jj,jl) = zindb* zs0ice(ji,jj,jl)353 a_i(ji,jj,jl) = zindb* zs0a (ji,jj,jl)354 e_s(ji,jj,1,jl) = zindb* zs0c0 (ji,jj,jl)344 v_s(ji,jj,jl) = rswitch * zs0sn (ji,jj,jl) 345 v_i(ji,jj,jl) = rswitch * zs0ice(ji,jj,jl) 346 a_i(ji,jj,jl) = rswitch * zs0a (ji,jj,jl) 347 e_s(ji,jj,1,jl) = rswitch * zs0c0 (ji,jj,jl) 355 348 ! Ice salinity and age 356 349 IF( num_sal == 2 ) THEN 357 350 smv_i(ji,jj,jl) = MAX( MIN( s_i_max * v_i(ji,jj,jl), zsmv ), s_i_min * v_i(ji,jj,jl) ) 358 351 ENDIF 359 oa_i(ji,jj,jl) = MAX( zindb* zs0oi(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ), 0._wp ) * a_i(ji,jj,jl)352 oa_i(ji,jj,jl) = MAX( rswitch * zs0oi(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ), 0._wp ) * a_i(ji,jj,jl) 360 353 361 354 ! Update fluxes … … 372 365 DO jj = 1, jpj 373 366 DO ji = 1, jpi 374 zindb= MAX( 0._wp , SIGN( 1._wp, zs0a(ji,jj,jl) - epsi10 ) )367 rswitch = MAX( 0._wp , SIGN( 1._wp, zs0a(ji,jj,jl) - epsi10 ) ) 375 368 zei = zs0e(ji,jj,jk,jl) 376 e_i(ji,jj,jk,jl) = zindb* MAX( 0._wp, zs0e(ji,jj,jk,jl) )369 e_i(ji,jj,jk,jl) = rswitch * MAX( 0._wp, zs0e(ji,jj,jk,jl) ) 377 370 ! Update fluxes 378 371 hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 … … 393 386 zsmv = smv_i(ji,jj,jl) 394 387 zes = e_s (ji,jj,1,jl) 395 zei = SUM( e_i(ji,jj, :,jl) )388 zei = SUM( e_i(ji,jj,1:nlay_i,jl) ) 396 389 zdv = v_i(ji,jj,jl) - zviold(ji,jj,jl) 397 390 !zda = a_i(ji,jj,jl) - zaiold(ji,jj,jl) 398 391 399 zindh = 1._wp392 rswitch = 1._wp 400 393 IF ( ( zdv > 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) .AND. SUM( zaiold(ji,jj,1:jpl) ) < 0.80 ) .OR. & 401 394 & ( zdv < 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) ) ) THEN 402 395 ht_i(ji,jj,jl) = MIN( zhimax(ji,jj,jl), hi_max(jl) ) 403 zindh =MAX( 0._wp, SIGN( 1._wp, ht_i(ji,jj,jl) - epsi20 ) )404 a_i(ji,jj,jl) = zindh * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi20 )396 rswitch = MAX( 0._wp, SIGN( 1._wp, ht_i(ji,jj,jl) - epsi20 ) ) 397 a_i(ji,jj,jl) = rswitch * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi20 ) 405 398 ELSE 406 399 ht_i(ji,jj,jl) = MAX( MIN( ht_i(ji,jj,jl), hi_max(jl) ), hi_max(jl-1) ) 407 zindh =MAX( 0._wp, SIGN( 1._wp, ht_i(ji,jj,jl) - epsi20 ) )408 a_i(ji,jj,jl) = zindh * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi20 )400 rswitch = MAX( 0._wp, SIGN( 1._wp, ht_i(ji,jj,jl) - epsi20 ) ) 401 a_i(ji,jj,jl) = rswitch * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi20 ) 409 402 ENDIF 410 403 411 ! small correction due to * zindh for a_i412 v_i (ji,jj,jl) = zindh * v_i (ji,jj,jl)413 v_s (ji,jj,jl) = zindh * v_s (ji,jj,jl)414 smv_i(ji,jj,jl) = zindh * smv_i(ji,jj,jl)415 e_s(ji,jj,1,jl) = zindh * e_s(ji,jj,1,jl)416 e_i(ji,jj, :,jl) = zindh * e_i(ji,jj,:,jl)404 ! small correction due to *rswitch for a_i 405 v_i (ji,jj,jl) = rswitch * v_i (ji,jj,jl) 406 v_s (ji,jj,jl) = rswitch * v_s (ji,jj,jl) 407 smv_i(ji,jj,jl) = rswitch * smv_i(ji,jj,jl) 408 e_s(ji,jj,1,jl) = rswitch * e_s(ji,jj,1,jl) 409 e_i(ji,jj,1:nlay_i,jl) = rswitch * e_i(ji,jj,1:nlay_i,jl) 417 410 418 411 ! Update mass fluxes … … 421 414 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice 422 415 hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 423 hfx_res(ji,jj) = hfx_res(ji,jj) + ( SUM( e_i(ji,jj,:,jl) ) - zei ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 424 416 hfx_res(ji,jj) = hfx_res(ji,jj) + ( SUM( e_i(ji,jj,1:nlay_i,jl) ) - zei ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 425 417 ENDIF 426 427 diag_trp_vi(ji,jj) = diag_trp_vi(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * r1_rdtice428 diag_trp_vs(ji,jj) = diag_trp_vs(ji,jj) + ( v_s(ji,jj,jl) - zvsold(ji,jj,jl) ) * r1_rdtice429 430 418 END DO 431 419 END DO … … 438 426 diag_trp_ei(ji,jj) = ( SUM( e_i(ji,jj,1:nlay_i,:) ) - zeiold(ji,jj) ) / area(ji,jj) * unit_fac * r1_rdtice 439 427 diag_trp_es(ji,jj) = ( SUM( e_s(ji,jj,1:nlay_s,:) ) - zesold(ji,jj) ) / area(ji,jj) * unit_fac * r1_rdtice 440 END DO 441 END DO 442 443 ! --- agglomerate variables (clem) ----------------- 428 429 diag_trp_vi(ji,jj) = SUM( v_i(ji,jj,:) - zviold(ji,jj,:) ) * r1_rdtice 430 diag_trp_vs(ji,jj) = SUM( v_s(ji,jj,:) - zvsold(ji,jj,:) ) * r1_rdtice 431 END DO 432 END DO 433 434 ! --- agglomerate variables ----------------- 444 435 vt_i (:,:) = 0._wp 445 436 vt_s (:,:) = 0._wp … … 462 453 DO ji = 1, jpi 463 454 ! open water = 1 if at_i=0 464 zindb= MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) )465 ato_i(ji,jj) = zindb + (1._wp - zindb) * zs0ow(ji,jj)455 rswitch = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) ) 456 ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * zs0ow(ji,jj) 466 457 END DO 467 458 END DO … … 506 497 CALL wrk_dealloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow, zeiold, zesold ) 507 498 CALL wrk_dealloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi ) 508 CALL wrk_dealloc( jpi, jpj, jkmax, jpl, zs0e )499 CALL wrk_dealloc( jpi, jpj, nlay_i+1, jpl, zs0e ) 509 500 510 501 CALL wrk_dealloc( jpi, jpj, jpl, zviold, zvsold, zaiold, zhimax ) ! clem
Note: See TracChangeset
for help on using the changeset viewer.
