- Timestamp:
- 2020-12-02T14:55:21+01:00 (3 years ago)
- Location:
- NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3
- Property svn:externals
-
old new 8 8 9 9 # SETTE 10 ^/utils/CI/sette@13 292sette10 ^/utils/CI/sette@13559 sette
-
- Property svn:externals
-
NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/ICE/icethd.F90
r13295 r13998 18 18 USE ice ! sea-ice: variables 19 19 !!gm list trop longue ==>>> why not passage en argument d'appel ? 20 USE sbc_oce , ONLY : sss_m, sst_m, e3t_m, utau, vtau, ssu_m, ssv_m, frq_m, qns_tot, qsr_tot,sprecip, ln_cpl20 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 22 & qml_ice, qcn_ice, qtr_ice_top … … 30 30 USE icethd_pnd ! sea-ice: melt ponds 31 31 USE iceitd ! sea-ice: remapping thickness distribution 32 USE icecor ! sea-ice: corrections 32 33 USE icetab ! sea-ice: 1D <==> 2D transformation 33 34 USE icevar ! sea-ice: operations … … 35 36 ! 36 37 USE in_out_manager ! I/O manager 38 USE iom ! I/O manager library 37 39 USE lib_mpp ! MPP library 38 40 USE lib_fortran ! fortran utilities (glob_sum + no signed zero) … … 51 53 LOGICAL :: ln_icedO ! activate ice growth in open-water (T) or not (F) 52 54 LOGICAL :: ln_icedS ! activate gravity drainage and flushing (T) or not (F) 55 LOGICAL :: ln_leadhfx ! heat in the leads is used to melt sea-ice before warming the ocean 56 57 !! for convergence tests 58 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztice_cvgerr, ztice_cvgstp 53 59 54 60 !! * Substitutions … … 86 92 ! 87 93 INTEGER :: ji, jj, jk, jl ! dummy loop indices 88 REAL(wp) :: zfric_u, zqld, zqfr, zqfr_neg 89 REAL(wp), PARAMETER :: zfric_umin = 0._wp 90 REAL(wp), PARAMETER :: zch = 0.0057_wp 91 REAL(wp), DIMENSION(jpi,jpj) :: zu_io, zv_io, zfric ! ice-ocean velocity (m/s) and frictional velocity (m2/s2)94 REAL(wp) :: zfric_u, zqld, zqfr, zqfr_neg, zqfr_pos 95 REAL(wp), PARAMETER :: zfric_umin = 0._wp ! lower bound for the friction velocity (cice value=5.e-04) 96 REAL(wp), PARAMETER :: zch = 0.0057_wp ! heat transfer coefficient 97 REAL(wp), DIMENSION(jpi,jpj) :: zu_io, zv_io, zfric, zvel ! ice-ocean velocity (m/s) and frictional velocity (m2/s2) 92 98 ! 93 99 !!------------------------------------------------------------------- … … 101 107 WRITE(numout,*) 'ice_thd: sea-ice thermodynamics' 102 108 WRITE(numout,*) '~~~~~~~' 109 ENDIF 110 111 ! convergence tests 112 IF( ln_zdf_chkcvg ) THEN 113 ALLOCATE( ztice_cvgerr(jpi,jpj,jpl) , ztice_cvgstp(jpi,jpj,jpl) ) 114 ztice_cvgerr = 0._wp ; ztice_cvgstp = 0._wp 103 115 ENDIF 104 116 … … 113 125 & ( zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj) & 114 126 & + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1) ) ) * tmask(ji,jj,1) 127 zvel(ji,jj) = 0.5_wp * SQRT( ( u_ice(ji-1,jj) + u_ice(ji,jj) ) * ( u_ice(ji-1,jj) + u_ice(ji,jj) ) + & 128 & ( v_ice(ji,jj-1) + v_ice(ji,jj) ) * ( v_ice(ji,jj-1) + v_ice(ji,jj) ) ) 115 129 END_2D 116 130 ELSE ! if no ice dynamics => transmit directly the atmospheric stress to the ocean … … 119 133 & ( utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj) & 120 134 & + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) ) * tmask(ji,jj,1) 135 zvel(ji,jj) = 0._wp 121 136 END_2D 122 137 ENDIF 123 CALL lbc_lnk ( 'icethd', zfric, 'T',1.0_wp )138 CALL lbc_lnk_multi( 'icethd', zfric, 'T', 1.0_wp, zvel, 'T', 1.0_wp ) 124 139 ! 125 140 !--------------------------------------------------------------------! … … 129 144 rswitch = tmask(ji,jj,1) * MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 0 if no ice 130 145 ! 131 ! ! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget132 ! ! practically no "direct lateral ablation"133 !134 ! ! net downward heat flux from the ice to the ocean, expressed as a function of ocean135 ! ! temperature and turbulent mixing (McPhee, 1992)136 !137 146 ! --- Energy received in the lead from atm-oce exchanges, zqld is defined everywhere (J.m-2) --- ! 138 147 zqld = tmask(ji,jj,1) * rDt_ice * & … … 140 149 & ( 1._wp - at_i_b(ji,jj) ) * qns_oce(ji,jj) + qemp_oce(ji,jj) ) 141 150 142 ! --- Energy needed to bring ocean surface layer until its freezing (mostly<0 but >0 if supercooling, J.m-2) --- ! 151 ! --- Energy needed to bring ocean surface layer until its freezing, zqfr is defined everywhere (J.m-2) --- ! 152 ! (mostly<0 but >0 if supercooling) 143 153 zqfr = rho0 * rcp * e3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) * tmask(ji,jj,1) ! both < 0 (t_bo < sst) and > 0 (t_bo > sst) 144 154 zqfr_neg = MIN( zqfr , 0._wp ) ! only < 0 145 146 ! --- Sensible ocean-to-ice heat flux (mostly>0 but <0 if supercooling, W/m2) 155 zqfr_pos = MAX( zqfr , 0._wp ) ! only > 0 156 157 ! --- Sensible ocean-to-ice heat flux (W/m2) --- ! 158 ! (mostly>0 but <0 if supercooling) 147 159 zfric_u = MAX( SQRT( zfric(ji,jj) ), zfric_umin ) 148 qsb_ice_bot(ji,jj) = rswitch * rho0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ! W.m-2 149 150 qsb_ice_bot(ji,jj) = rswitch * MIN( qsb_ice_bot(ji,jj), - zqfr_neg * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) ) 160 qsb_ice_bot(ji,jj) = rswitch * rho0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) 161 151 162 ! upper bound for qsb_ice_bot: the heat retrieved from the ocean must be smaller than the heat necessary to reach 152 163 ! the freezing point, so that we do not have SST < T_freeze 153 ! This implies: - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rtdice ) - zqfr >= 0 154 155 !-- Energy Budget of the leads (J.m-2), source of ice growth in open water. Must be < 0 to form ice 156 qlead(ji,jj) = MIN( 0._wp , zqld - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rDt_ice ) - zqfr ) 157 158 ! If there is ice and leads are warming => transfer energy from the lead budget and use it for bottom melting 159 ! If the grid cell is fully covered by ice (no leads) => transfer energy from the lead budget to the ice bottom budget 160 IF( ( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) .OR. at_i(ji,jj) >= (1._wp - epsi10) ) THEN 161 fhld (ji,jj) = rswitch * zqld * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90 164 ! This implies: qsb_ice_bot(ji,jj) * at_i(ji,jj) * rtdice <= - zqfr_neg 165 ! The following formulation is ok for both normal conditions and supercooling 166 qsb_ice_bot(ji,jj) = rswitch * MIN( qsb_ice_bot(ji,jj), - zqfr_neg * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) ) 167 168 ! --- Energy Budget of the leads (qlead, J.m-2) --- ! 169 ! qlead is the energy received from the atm. in the leads. 170 ! If warming (zqld >= 0), then the energy in the leads is used to melt ice (bottom melting) => fhld (W/m2) 171 ! If cooling (zqld < 0), then the energy in the leads is used to grow ice in open water => qlead (J.m-2) 172 IF( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) THEN 173 ! upper bound for fhld: fhld should be equal to zqld 174 ! but we have to make sure that this heat will not make the sst drop below the freezing point 175 ! so the max heat that can be pulled out of the ocean is zqld - qsb - zqfr_pos 176 ! The following formulation is ok for both normal conditions and supercooling 177 fhld (ji,jj) = rswitch * MAX( 0._wp, ( zqld - zqfr_pos ) * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) & ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90 178 & - qsb_ice_bot(ji,jj) ) 162 179 qlead(ji,jj) = 0._wp 163 180 ELSE 164 181 fhld (ji,jj) = 0._wp 182 ! upper bound for qlead: qlead should be equal to zqld 183 ! but before using this heat for ice formation, we suppose that the ocean cools down till the freezing point. 184 ! The energy for this cooling down is zqfr. Also some heat will be removed from the ocean from turbulent fluxes (qsb) 185 ! and freezing point is reached if zqfr = zqld - qsb*a/dt 186 ! so the max heat that can be pulled out of the ocean is zqld - qsb - zqfr 187 ! The following formulation is ok for both normal conditions and supercooling 188 qlead(ji,jj) = MIN( 0._wp , zqld - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rDt_ice ) - zqfr ) 165 189 ENDIF 166 190 ! 167 ! Net heat flux on top of the ice-ocean [W.m-2] 168 ! --------------------------------------------- 169 qt_atm_oi(ji,jj) = qns_tot(ji,jj) + qsr_tot(ji,jj) 191 ! If ice is landfast and ice concentration reaches its max 192 ! => stop ice formation in open water 193 IF( zvel(ji,jj) <= 5.e-04_wp .AND. at_i(ji,jj) >= rn_amax_2d(ji,jj)-epsi06 ) qlead(ji,jj) = 0._wp 194 ! 195 ! If the grid cell is almost fully covered by ice (no leads) 196 ! => stop ice formation in open water 197 IF( at_i(ji,jj) >= (1._wp - epsi10) ) qlead(ji,jj) = 0._wp 198 ! 199 ! If ln_leadhfx is false 200 ! => do not use energy of the leads to melt sea-ice 201 IF( .NOT.ln_leadhfx ) fhld(ji,jj) = 0._wp 202 ! 170 203 END_2D 171 204 … … 178 211 ENDIF 179 212 180 ! ---------------------------------------------------------------------181 ! Net heat flux on top of the ocean after ice thermo (1st step) [W.m-2]182 ! ---------------------------------------------------------------------183 ! First step here : non solar + precip - qlead - qsensible184 ! Second step in icethd_dh : heat remaining if total melt (zq_rema)185 ! Third step in iceupdate.F90 : heat from ice-ocean mass exchange (zf_mass) + solar186 qt_oce_ai(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + qemp_oce(:,:) & ! Non solar heat flux received by the ocean187 & - qlead(:,:) * r1_Dt_ice & ! heat flux taken from the ocean where there is open water ice formation188 & - at_i (:,:) * qsb_ice_bot(:,:) & ! heat flux taken by sensible flux189 & - at_i (:,:) * fhld (:,:) ! heat flux taken during bottom growth/melt190 ! ! (fhld should be 0 while bott growth)191 213 !-------------------------------------------------------------------------------------------! 192 214 ! Thermodynamic computation (only on grid points covered by ice) => loop over ice categories … … 208 230 ! ! --- & Change units of e_i, e_s from J/m2 to J/m3 --- ! 209 231 ! 210 s_i_new (1:npti) = 0._wp ; dh_s_tot(1:npti) = 0._wp ! --- some init --- ! (important to have them here)232 s_i_new (1:npti) = 0._wp ; dh_s_tot(1:npti) = 0._wp ! --- some init --- ! (important to have them here) 211 233 dh_i_sum (1:npti) = 0._wp ; dh_i_bom(1:npti) = 0._wp ; dh_i_itm (1:npti) = 0._wp 212 234 dh_i_sub (1:npti) = 0._wp ; dh_i_bog(1:npti) = 0._wp … … 218 240 CALL ice_thd_dh ! Ice-Snow thickness 219 241 CALL ice_thd_pnd ! Melt ponds formation 220 CALL ice_thd_ent( e_i_1d(1:npti,:) , .true.) ! Ice enthalpy remapping242 CALL ice_thd_ent( e_i_1d(1:npti,:) ) ! Ice enthalpy remapping 221 243 ENDIF 222 244 CALL ice_thd_sal( ln_icedS ) ! --- Ice salinity --- ! … … 241 263 ! 242 264 IF( ln_icedO ) CALL ice_thd_do ! --- Frazil ice growth in leads --- ! 265 ! 266 CALL ice_cor( kt , 2 ) ! --- Corrections --- ! 267 ! 268 oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rdt_ice ! ice natural aging incrementation 269 ! 270 ! convergence tests 271 IF( ln_zdf_chkcvg ) THEN 272 CALL iom_put( 'tice_cvgerr', ztice_cvgerr ) ; DEALLOCATE( ztice_cvgerr ) 273 CALL iom_put( 'tice_cvgstp', ztice_cvgstp ) ; DEALLOCATE( ztice_cvgstp ) 274 ENDIF 243 275 ! 244 276 ! controls … … 347 379 CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip (:,:,kl) ) 348 380 CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip (:,:,kl) ) 349 CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_frac_1d(1:npti), a_ip_frac(:,:,kl) )381 CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d (1:npti), h_il (:,:,kl) ) 350 382 ! 351 383 CALL tab_2d_1d( npti, nptidx(1:npti), qprec_ice_1d (1:npti), qprec_ice ) … … 399 431 CALL tab_2d_1d( npti, nptidx(1:npti), hfx_res_1d (1:npti), hfx_res ) 400 432 CALL tab_2d_1d( npti, nptidx(1:npti), hfx_err_dif_1d(1:npti), hfx_err_dif ) 401 CALL tab_2d_1d( npti, nptidx(1:npti), hfx_err_rem_1d(1:npti), hfx_err_rem )402 CALL tab_2d_1d( npti, nptidx(1:npti), qt_oce_ai_1d (1:npti), qt_oce_ai )403 433 ! 404 434 ! ocean surface fields 405 435 CALL tab_2d_1d( npti, nptidx(1:npti), sst_1d(1:npti), sst_m ) 406 436 CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d(1:npti), sss_m ) 437 CALL tab_2d_1d( npti, nptidx(1:npti), frq_m_1d(1:npti), frq_m ) 407 438 ! 408 439 ! to update ice age … … 434 465 sv_i_1d(1:npti) = s_i_1d (1:npti) * v_i_1d (1:npti) 435 466 v_ip_1d(1:npti) = h_ip_1d(1:npti) * a_ip_1d(1:npti) 467 v_il_1d(1:npti) = h_il_1d(1:npti) * a_ip_1d(1:npti) 436 468 oa_i_1d(1:npti) = o_i_1d (1:npti) * a_i_1d (1:npti) 437 469 … … 453 485 CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip (:,:,kl) ) 454 486 CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip (:,:,kl) ) 455 CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_frac_1d(1:npti), a_ip_frac(:,:,kl) )487 CALL tab_1d_2d( npti, nptidx(1:npti), h_il_1d (1:npti), h_il (:,:,kl) ) 456 488 ! 457 489 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni ) … … 491 523 CALL tab_1d_2d( npti, nptidx(1:npti), hfx_res_1d (1:npti), hfx_res ) 492 524 CALL tab_1d_2d( npti, nptidx(1:npti), hfx_err_dif_1d(1:npti), hfx_err_dif ) 493 CALL tab_1d_2d( npti, nptidx(1:npti), hfx_err_rem_1d(1:npti), hfx_err_rem )494 CALL tab_1d_2d( npti, nptidx(1:npti), qt_oce_ai_1d (1:npti), qt_oce_ai )495 525 ! 496 526 CALL tab_1d_2d( npti, nptidx(1:npti), qns_ice_1d (1:npti), qns_ice (:,:,kl) ) … … 508 538 CALL tab_1d_2d( npti, nptidx(1:npti), sv_i_1d(1:npti), sv_i(:,:,kl) ) 509 539 CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d(1:npti), v_ip(:,:,kl) ) 540 CALL tab_1d_2d( npti, nptidx(1:npti), v_il_1d(1:npti), v_il(:,:,kl) ) 510 541 CALL tab_1d_2d( npti, nptidx(1:npti), oa_i_1d(1:npti), oa_i(:,:,kl) ) 542 ! check convergence of heat diffusion scheme 543 IF( ln_zdf_chkcvg ) THEN 544 CALL tab_1d_2d( npti, nptidx(1:npti), tice_cvgerr_1d(1:npti), ztice_cvgerr(:,:,kl) ) 545 CALL tab_1d_2d( npti, nptidx(1:npti), tice_cvgstp_1d(1:npti), ztice_cvgstp(:,:,kl) ) 546 ENDIF 511 547 ! 512 548 END SELECT … … 529 565 INTEGER :: ios ! Local integer output status for namelist read 530 566 !! 531 NAMELIST/namthd/ ln_icedH, ln_icedA, ln_icedO, ln_icedS 567 NAMELIST/namthd/ ln_icedH, ln_icedA, ln_icedO, ln_icedS, ln_leadhfx 532 568 !!------------------------------------------------------------------- 533 569 ! … … 543 579 WRITE(numout,*) '~~~~~~~~~~~~' 544 580 WRITE(numout,*) ' Namelist namthd:' 545 WRITE(numout,*) ' activate ice thick change from top/bot (T) or not (F) ln_icedH = ', ln_icedH 546 WRITE(numout,*) ' activate lateral melting (T) or not (F) ln_icedA = ', ln_icedA 547 WRITE(numout,*) ' activate ice growth in open-water (T) or not (F) ln_icedO = ', ln_icedO 548 WRITE(numout,*) ' activate gravity drainage and flushing (T) or not (F) ln_icedS = ', ln_icedS 581 WRITE(numout,*) ' activate ice thick change from top/bot (T) or not (F) ln_icedH = ', ln_icedH 582 WRITE(numout,*) ' activate lateral melting (T) or not (F) ln_icedA = ', ln_icedA 583 WRITE(numout,*) ' activate ice growth in open-water (T) or not (F) ln_icedO = ', ln_icedO 584 WRITE(numout,*) ' activate gravity drainage and flushing (T) or not (F) ln_icedS = ', ln_icedS 585 WRITE(numout,*) ' heat in the leads is used to melt sea-ice before warming the ocean ln_leadhfx = ', ln_leadhfx 549 586 ENDIF 550 587 !
Note: See TracChangeset
for help on using the changeset viewer.