- Timestamp:
- 2015-07-21T13:25:36+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
r5146 r5621 22 22 USE phycst ! physical constants 23 23 USE dom_oce ! ocean space and time domain variables 24 USE oce , ONLY : fraqsr_1lev25 24 USE ice ! LIM: sea-ice variables 26 25 USE sbc_oce ! Surface boundary condition: ocean fields … … 28 27 USE thd_ice ! LIM thermodynamic sea-ice variables 29 28 USE dom_ice ! LIM sea-ice domain 30 USE domvvl ! domain: variable volume level31 29 USE limthd_dif ! LIM: thermodynamics, vertical diffusion 32 30 USE limthd_dh ! LIM: thermodynamics, ice and snow thickness variation … … 50 48 PRIVATE 51 49 52 PUBLIC lim_thd ! called by limstp module53 PUBLIC lim_thd_init ! called by sbc_lim_init50 PUBLIC lim_thd ! called by limstp module 51 PUBLIC lim_thd_init ! called by sbc_lim_init 54 52 55 53 !! * Substitutions … … 89 87 REAL(wp) :: zfric_u, zqld, zqfr 90 88 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 91 REAL(wp), PARAMETER :: zfric_umin = 0._wp ! lower bound for the friction velocity (cice value=5.e-04) 92 REAL(wp), PARAMETER :: zch = 0.0057_wp ! heat transfer coefficient 93 ! 94 REAL(wp), POINTER, DIMENSION(:,:) :: zqsr, zqns 89 REAL(wp), PARAMETER :: zfric_umin = 0._wp ! lower bound for the friction velocity (cice value=5.e-04) 90 REAL(wp), PARAMETER :: zch = 0.0057_wp ! heat transfer coefficient 91 ! 95 92 !!------------------------------------------------------------------- 96 CALL wrk_alloc( jpi, jpj, zqsr, zqns )97 93 98 94 IF( nn_timing == 1 ) CALL timing_start('limthd') … … 101 97 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 102 98 99 CALL lim_var_glo2eqv 103 100 !------------------------------------------------------------------------! 104 101 ! 1) Initialization of some variables ! … … 135 132 ! 2) Partial computation of forcing for the thermodynamic sea ice model. ! 136 133 !-----------------------------------------------------------------------------! 137 138 !--- Ocean solar and non solar fluxes to be used in zqld139 IF ( .NOT. lk_cpl ) THEN ! --- forced case, fluxes to the lead are the same as over the ocean140 !141 zqsr(:,:) = qsr(:,:) ; zqns(:,:) = qns(:,:)142 !143 ELSE ! --- coupled case, fluxes to the lead are total - intercepted144 !145 zqsr(:,:) = qsr_tot(:,:) ; zqns(:,:) = qns_tot(:,:)146 !147 DO jl = 1, jpl148 DO jj = 1, jpj149 DO ji = 1, jpi150 zqsr(ji,jj) = zqsr(ji,jj) - qsr_ice(ji,jj,jl) * a_i_b(ji,jj,jl)151 zqns(ji,jj) = zqns(ji,jj) - qns_ice(ji,jj,jl) * a_i_b(ji,jj,jl)152 END DO153 END DO154 END DO155 !156 ENDIF157 158 134 DO jj = 1, jpj 159 135 DO ji = 1, jpi … … 166 142 ! ! temperature and turbulent mixing (McPhee, 1992) 167 143 ! 168 169 144 ! --- Energy received in the lead, zqld is defined everywhere (J.m-2) --- ! 170 ! REMARK valid at least in forced mode from clem 171 ! precip is included in qns but not in qns_ice 172 IF ( lk_cpl ) THEN 173 zqld = tmask(ji,jj,1) * rdt_ice * & 174 & ( zqsr(ji,jj) * fraqsr_1lev(ji,jj) + zqns(ji,jj) & ! pfrld already included in coupled mode 175 & + ( pfrld(ji,jj)**rn_betas - pfrld(ji,jj) ) * sprecip(ji,jj) * & ! heat content of precip 176 & ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rt0 ) - lfus ) & 177 & + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rt0 ) ) 178 ELSE 179 zqld = tmask(ji,jj,1) * rdt_ice * & 180 & ( pfrld(ji,jj) * ( zqsr(ji,jj) * fraqsr_1lev(ji,jj) + zqns(ji,jj) ) & 181 & + ( pfrld(ji,jj)**rn_betas - pfrld(ji,jj) ) * sprecip(ji,jj) * & ! heat content of precip 182 & ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rt0 ) - lfus ) & 183 & + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rt0 ) ) 184 ENDIF 145 zqld = tmask(ji,jj,1) * rdt_ice * & 146 & ( pfrld(ji,jj) * qsr_oce(ji,jj) * frq_m(ji,jj) + pfrld(ji,jj) * qns_oce(ji,jj) + qemp_oce(ji,jj) ) 185 147 186 148 ! --- Energy needed to bring ocean surface layer until its freezing (<0, J.m-2) --- ! … … 209 171 ! Net heat flux on top of ice-ocean [W.m-2] 210 172 ! ----------------------------------------- 211 ! First step here : heat flux at the ocean surface + precip 212 ! Second step below : heat flux at the ice surface (after limthd_dif) 213 hfx_in(ji,jj) = hfx_in(ji,jj) & 214 ! heat flux above the ocean 215 & + pfrld(ji,jj) * ( zqns(ji,jj) + zqsr(ji,jj) ) & 216 ! latent heat of precip (note that precip is included in qns but not in qns_ice) 217 & + ( 1._wp - pfrld(ji,jj) ) * sprecip(ji,jj) * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rt0 ) - lfus ) & 218 & + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rt0 ) 173 hfx_in(ji,jj) = qns_tot(ji,jj) + qsr_tot(ji,jj) 219 174 220 175 ! ----------------------------------------------------------------------------- 221 ! Net heat flux that is retroceded to the ocean or taken from the ocean[W.m-2]176 ! Net heat flux on top of the ocean after ice thermo (1st step) [W.m-2] 222 177 ! ----------------------------------------------------------------------------- 223 178 ! First step here : non solar + precip - qlead - qturb 224 179 ! Second step in limthd_dh : heat remaining if total melt (zq_rema) 225 180 ! Third step in limsbc : heat from ice-ocean mass exchange (zf_mass) + solar 226 hfx_out(ji,jj) = hfx_out(ji,jj) & 227 ! Non solar heat flux received by the ocean 228 & + pfrld(ji,jj) * qns(ji,jj) & 229 ! latent heat of precip (note that precip is included in qns but not in qns_ice) 230 & + ( pfrld(ji,jj)**rn_betas - pfrld(ji,jj) ) * sprecip(ji,jj) & 231 & * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rt0 ) - lfus ) & 232 & + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rt0 ) & 233 ! heat flux taken from the ocean where there is open water ice formation 234 & - qlead(ji,jj) * r1_rdtice & 235 ! heat flux taken from the ocean during bottom growth/melt (fhld should be 0 while bott growth) 236 & - at_i(ji,jj) * fhtur(ji,jj) & 237 & - at_i(ji,jj) * fhld(ji,jj) 238 181 hfx_out(ji,jj) = pfrld(ji,jj) * qns_oce(ji,jj) + qemp_oce(ji,jj) & ! Non solar heat flux received by the ocean 182 & - qlead(ji,jj) * r1_rdtice & ! heat flux taken from the ocean where there is open water ice formation 183 & - at_i(ji,jj) * fhtur(ji,jj) & ! heat flux taken by turbulence 184 & - at_i(ji,jj) * fhld(ji,jj) ! heat flux taken during bottom growth/melt 185 ! (fhld should be 0 while bott growth) 239 186 END DO 240 187 END DO … … 311 258 ! --- lateral melting if monocat --- ! 312 259 !------------------------------------! 313 IF ( ( ( nn_monocat == 1 ) .OR. ( nn_monocat == 4 ) ) .AND. ( jpl == 1 )) THEN260 IF ( ( nn_monocat == 1 .OR. nn_monocat == 4 ) .AND. jpl == 1 ) THEN 314 261 CALL lim_thd_lam( 1, nbpb ) 315 262 END IF … … 324 271 ENDIF 325 272 ! 326 END DO 273 END DO !jl 327 274 328 275 !------------------------------------------------------------------------------! … … 350 297 END DO 351 298 352 !------------------------353 ! Ice natural aging354 !------------------------355 oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rdt_ice /rday356 357 299 !---------------------------------- 358 300 ! Change thickness to volume 359 301 !---------------------------------- 360 CALL lim_var_eqv2glo 302 v_i(:,:,:) = ht_i(:,:,:) * a_i(:,:,:) 303 v_s(:,:,:) = ht_s(:,:,:) * a_i(:,:,:) 304 smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:) 305 306 ! update ice age (in case a_i changed, i.e. becomes 0 or lateral melting in monocat) 307 DO jl = 1, jpl 308 DO jj = 1, jpj 309 DO ji = 1, jpi 310 rswitch = MAX( 0._wp , SIGN( 1._wp, a_i_b(ji,jj,jl) - epsi10 ) ) 311 oa_i(ji,jj,jl) = rswitch * oa_i(ji,jj,jl) * a_i(ji,jj,jl) / MAX( a_i_b(ji,jj,jl), epsi10 ) 312 END DO 313 END DO 314 END DO 361 315 362 316 CALL lim_var_zapsmall 317 363 318 !-------------------------------------------- 364 319 ! Diagnostic thermodynamic growth rates … … 399 354 ! 400 355 ! 401 CALL wrk_dealloc( jpi, jpj, zqsr, zqns )402 403 356 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 357 404 358 !------------------------------------------------------------------------------| 405 359 ! 6) Transport of ice between thickness categories. | 406 360 !------------------------------------------------------------------------------| 361 ! Given thermodynamic growth rates, transport ice between thickness categories. 407 362 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 408 363 409 ! Given thermodynamic growth rates, transport ice between thickness categories. 410 IF( jpl > 1 ) CALL lim_itd_th_rem( 1, jpl, kt ) 411 ! 412 CALL lim_var_glo2eqv ! only for info 413 CALL lim_var_agg(1) 364 IF( jpl > 1 ) CALL lim_itd_th_rem( 1, jpl, kt ) 414 365 415 366 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 367 416 368 !------------------------------------------------------------------------------| 417 369 ! 7) Add frazil ice growing in leads. 418 370 !------------------------------------------------------------------------------| 419 371 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 372 420 373 CALL lim_thd_lac 421 CALL lim_var_glo2eqv ! only for info422 374 423 ! conservation test424 375 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 425 376 426 IF(ln_ctl) THEN ! Control print 377 ! Control print 378 IF(ln_ctl) THEN 379 CALL lim_var_glo2eqv 380 427 381 CALL prt_ctl_info(' ') 428 382 CALL prt_ctl_info(' - Cell values : ') … … 460 414 END SUBROUTINE lim_thd 461 415 416 462 417 SUBROUTINE lim_thd_temp( kideb, kiut ) 463 418 !!----------------------------------------------------------------------- … … 503 458 REAL(wp) :: zhi_bef ! ice thickness before thermo 504 459 REAL(wp) :: zdh_mel, zda_mel ! net melting 505 REAL(wp) :: zv ! ice volume460 REAL(wp) :: zvi, zvs ! ice/snow volumes 506 461 507 462 DO ji = kideb, kiut 508 463 zdh_mel = MIN( 0._wp, dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) ) 509 IF( zdh_mel < 0._wp ) THEN 510 zv = a_i_1d(ji) * ht_i_1d(ji) 464 IF( zdh_mel < 0._wp .AND. a_i_1d(ji) > 0._wp ) THEN 465 zvi = a_i_1d(ji) * ht_i_1d(ji) 466 zvs = a_i_1d(ji) * ht_s_1d(ji) 511 467 ! lateral melting = concentration change 512 468 zhi_bef = ht_i_1d(ji) - zdh_mel 513 zda_mel = a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi10 ) ) 514 a_i_1d(ji) = MAX( 0._wp, a_i_1d(ji) + zda_mel ) 515 ! adjust thickness 516 rswitch = MAX( 0._wp , SIGN( 1._wp , a_i_1d(ji) - epsi20 ) ) 517 ht_i_1d(ji) = rswitch * zv / MAX( a_i_1d(ji), epsi20 ) 469 rswitch = MAX( 0._wp , SIGN( 1._wp , zhi_bef - epsi20 ) ) 470 zda_mel = rswitch * a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi20 ) ) 471 a_i_1d(ji) = MAX( epsi20, a_i_1d(ji) + zda_mel ) 472 ! adjust thickness 473 ht_i_1d(ji) = zvi / a_i_1d(ji) 474 ht_s_1d(ji) = zvs / a_i_1d(ji) 518 475 ! retrieve total concentration 519 476 at_i_1d(ji) = a_i_1d(ji) … … 556 513 END DO 557 514 558 CALL tab_2d_1d( nbpb, tatm_ice_1d(1:nbpb), tatm_ice(:,:), jpi, jpj, npb(1:nbpb) )515 CALL tab_2d_1d( nbpb, qprec_ice_1d(1:nbpb), qprec_ice(:,:) , jpi, jpj, npb(1:nbpb) ) 559 516 CALL tab_2d_1d( nbpb, qsr_ice_1d (1:nbpb), qsr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 560 517 CALL tab_2d_1d( nbpb, fr1_i0_1d (1:nbpb), fr1_i0 , jpi, jpj, npb(1:nbpb) ) … … 562 519 CALL tab_2d_1d( nbpb, qns_ice_1d (1:nbpb), qns_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 563 520 CALL tab_2d_1d( nbpb, ftr_ice_1d (1:nbpb), ftr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 564 IF( .NOT. lk_cpl ) THEN 565 CALL tab_2d_1d( nbpb, qla_ice_1d (1:nbpb), qla_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 566 CALL tab_2d_1d( nbpb, dqla_ice_1d(1:nbpb), dqla_ice(:,:,jl), jpi, jpj, npb(1:nbpb) ) 567 ENDIF 521 CALL tab_2d_1d( nbpb, evap_ice_1d (1:nbpb), evap_ice(:,:,jl), jpi, jpj, npb(1:nbpb) ) 568 522 CALL tab_2d_1d( nbpb, dqns_ice_1d(1:nbpb), dqns_ice(:,:,jl), jpi, jpj, npb(1:nbpb) ) 569 523 CALL tab_2d_1d( nbpb, t_bo_1d (1:nbpb), t_bo , jpi, jpj, npb(1:nbpb) ) … … 656 610 CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qns_ice_1d(1:nbpb) , jpi, jpj) 657 611 CALL tab_1d_2d( nbpb, ftr_ice(:,:,jl), npb, ftr_ice_1d(1:nbpb) , jpi, jpj ) 658 612 ! 659 613 END SELECT 660 614 … … 676 630 INTEGER :: ios ! Local integer output status for namelist read 677 631 NAMELIST/namicethd/ rn_hnewice, ln_frazil, rn_maxfrazb, rn_vfrazb, rn_Cfrazb, & 678 & rn_himin, parsub,rn_betas, rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon, &632 & rn_himin, rn_betas, rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon, & 679 633 & nn_monocat, ln_it_qnsice 680 634 !!------------------------------------------------------------------- … … 700 654 ENDIF 701 655 702 IF( lk_cpl .AND. parsub /= 0.0 ) CALL ctl_stop( 'In coupled mode, use parsub = 0. or send dqla' )703 656 ! 704 657 IF(lwp) THEN ! control print … … 712 665 WRITE(numout,*)' minimum ice thickness rn_himin = ', rn_himin 713 666 WRITE(numout,*)' numerical carac. of the scheme for diffusion in ice ' 714 WRITE(numout,*)' switch for snow sublimation (=1) or not (=0) parsub = ', parsub715 667 WRITE(numout,*)' coefficient for ice-lead partition of snowfall rn_betas = ', rn_betas 716 668 WRITE(numout,*)' extinction radiation parameter in sea ice rn_kappa_i = ', rn_kappa_i
Note: See TracChangeset
for help on using the changeset viewer.
