- Timestamp:
- 2021-12-03T20:32:50+01:00 (3 years ago)
- Location:
- NEMO/branches/2021/dev_r14318_RK3_stage1
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/ICE/icethd.F90 (modified) (8 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r14318_RK3_stage1
- Property svn:externals
-
old new 9 9 10 10 # SETTE 11 ^/utils/CI/sette@14244 sette 11 ^/utils/CI/sette@HEAD sette 12
-
- Property svn:externals
-
NEMO/branches/2021/dev_r14318_RK3_stage1/src/ICE/icethd.F90
r14072 r15574 48 48 PUBLIC ice_thd_init ! called by ice_init 49 49 50 !!** namelist (namthd) **51 LOGICAL :: ln_icedH ! activate ice thickness change from growing/melting (T) or not (F)52 LOGICAL :: ln_icedA ! activate lateral melting param. (T) or not (F)53 LOGICAL :: ln_icedO ! activate ice growth in open-water (T) or not (F)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 ocean56 57 50 !! for convergence tests 58 51 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztice_cvgerr, ztice_cvgstp … … 92 85 ! 93 86 INTEGER :: ji, jj, jk, jl ! dummy loop indices 94 REAL(wp) :: zfric_u, zqld, zqfr, zqfr_neg, zqfr_pos95 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 coefficient97 REAL(wp), DIMENSION(jpi,jpj) :: zu_io, zv_io, zfric, zvel ! ice-ocean velocity (m/s) and frictional velocity (m2/s2)98 !99 87 !!------------------------------------------------------------------- 100 88 ! controls … … 114 102 ztice_cvgerr = 0._wp ; ztice_cvgstp = 0._wp 115 103 ENDIF 116 117 !---------------------------------------------! 118 ! computation of friction velocity at T points 119 !---------------------------------------------! 120 IF( ln_icedyn ) THEN 121 zu_io(:,:) = u_ice(:,:) - ssu_m(:,:) 122 zv_io(:,:) = v_ice(:,:) - ssv_m(:,:) 123 DO_2D( 0, 0, 0, 0 ) 124 zfric(ji,jj) = rn_cio * ( 0.5_wp * & 125 & ( zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj) & 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) ) ) 129 END_2D 130 ELSE ! if no ice dynamics => transmit directly the atmospheric stress to the ocean 131 DO_2D( 0, 0, 0, 0 ) 132 zfric(ji,jj) = r1_rho0 * SQRT( 0.5_wp * & 133 & ( utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj) & 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 136 END_2D 137 ENDIF 138 CALL lbc_lnk_multi( 'icethd', zfric, 'T', 1.0_wp, zvel, 'T', 1.0_wp ) 139 ! 140 !--------------------------------------------------------------------! 141 ! Partial computation of forcing for the thermodynamic sea ice model 142 !--------------------------------------------------------------------! 143 DO_2D( 1, 1, 1, 1 ) 144 rswitch = tmask(ji,jj,1) * MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 0 if no ice 145 ! 146 ! --- Energy received in the lead from atm-oce exchanges, zqld is defined everywhere (J.m-2) --- ! 147 zqld = tmask(ji,jj,1) * rDt_ice * & 148 & ( ( 1._wp - at_i_b(ji,jj) ) * qsr_oce(ji,jj) * frq_m(ji,jj) + & 149 & ( 1._wp - at_i_b(ji,jj) ) * qns_oce(ji,jj) + qemp_oce(ji,jj) ) 150 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) 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) 154 zqfr_neg = MIN( zqfr , 0._wp ) ! only < 0 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) 159 zfric_u = MAX( SQRT( zfric(ji,jj) ), zfric_umin ) 160 qsb_ice_bot(ji,jj) = rswitch * rho0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) 161 162 ! upper bound for qsb_ice_bot: the heat retrieved from the ocean must be smaller than the heat necessary to reach 163 ! the freezing point, so that we do not have SST < T_freeze 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 ! If conditions are always supercooled (such as at the mouth of ice-shelves), then ice grows continuously 169 ! ==> stop ice formation by artificially setting up the turbulent fluxes to 0 when volume > 20m (arbitrary) 170 IF( ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) > 0._wp .AND. vt_i(ji,jj) >= 20._wp ) THEN 171 zqfr = 0._wp 172 zqfr_pos = 0._wp 173 qsb_ice_bot(ji,jj) = 0._wp 174 ENDIF 175 ! 176 ! --- Energy Budget of the leads (qlead, J.m-2) --- ! 177 ! qlead is the energy received from the atm. in the leads. 178 ! If warming (zqld >= 0), then the energy in the leads is used to melt ice (bottom melting) => fhld (W/m2) 179 ! If cooling (zqld < 0), then the energy in the leads is used to grow ice in open water => qlead (J.m-2) 180 IF( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) THEN 181 ! upper bound for fhld: fhld should be equal to zqld 182 ! but we have to make sure that this heat will not make the sst drop below the freezing point 183 ! so the max heat that can be pulled out of the ocean is zqld - qsb - zqfr_pos 184 ! The following formulation is ok for both normal conditions and supercooling 185 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 186 & - qsb_ice_bot(ji,jj) ) 187 qlead(ji,jj) = 0._wp 188 ELSE 189 fhld (ji,jj) = 0._wp 190 ! upper bound for qlead: qlead should be equal to zqld 191 ! but before using this heat for ice formation, we suppose that the ocean cools down till the freezing point. 192 ! The energy for this cooling down is zqfr. Also some heat will be removed from the ocean from turbulent fluxes (qsb) 193 ! and freezing point is reached if zqfr = zqld - qsb*a/dt 194 ! so the max heat that can be pulled out of the ocean is zqld - qsb - zqfr 195 ! The following formulation is ok for both normal conditions and supercooling 196 qlead(ji,jj) = MIN( 0._wp , zqld - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rDt_ice ) - zqfr ) 197 ENDIF 198 ! 199 ! If ice is landfast and ice concentration reaches its max 200 ! => stop ice formation in open water 201 IF( zvel(ji,jj) <= 5.e-04_wp .AND. at_i(ji,jj) >= rn_amax_2d(ji,jj)-epsi06 ) qlead(ji,jj) = 0._wp 202 ! 203 ! If the grid cell is almost fully covered by ice (no leads) 204 ! => stop ice formation in open water 205 IF( at_i(ji,jj) >= (1._wp - epsi10) ) qlead(ji,jj) = 0._wp 206 ! 207 ! If ln_leadhfx is false 208 ! => do not use energy of the leads to melt sea-ice 209 IF( .NOT.ln_leadhfx ) fhld(ji,jj) = 0._wp 210 ! 211 END_2D 212 213 ! In case we bypass open-water ice formation 214 IF( .NOT. ln_icedO ) qlead(:,:) = 0._wp 215 ! In case we bypass growing/melting from top and bottom 216 IF( .NOT. ln_icedH ) THEN 217 qsb_ice_bot(:,:) = 0._wp 218 fhld (:,:) = 0._wp 219 ENDIF 220 104 ! 105 CALL ice_thd_frazil !--- frazil ice: collection thickness (ht_i_new) & fraction of frazil (fraz_frac) 106 ! 221 107 !-------------------------------------------------------------------------------------------! 222 108 ! Thermodynamic computation (only on grid points covered by ice) => loop over ice categories … … 226 112 ! select ice covered grid points 227 113 npti = 0 ; nptidx(:) = 0 228 DO_2D( 1, 1, 1, 1)114 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 229 115 IF ( a_i(ji,jj,jl) > epsi10 ) THEN 230 116 npti = npti + 1 … … 268 154 ! 269 155 IF ( ln_pnd .AND. ln_icedH ) & 270 & CALL ice_thd_pnd ! --- Melt ponds 156 & CALL ice_thd_pnd ! --- Melt ponds --- ! 271 157 ! 272 158 IF( jpl > 1 ) CALL ice_itd_rem( kt ) ! --- Transport ice between thickness categories --- ! … … 276 162 CALL ice_cor( kt , 2 ) ! --- Corrections --- ! 277 163 ! 278 oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rDt_ice ! ice natural aging incrementation 164 oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rDt_ice ! --- Ice natural aging incrementation 165 ! 166 DO_2D( 0, 0, 0, 0 ) ! --- Ice velocity corrections 167 IF( at_i(ji,jj) == 0._wp ) THEN ! if ice has melted 168 IF( at_i(ji+1,jj) == 0._wp ) u_ice(ji ,jj) = 0._wp ! right side 169 IF( at_i(ji-1,jj) == 0._wp ) u_ice(ji-1,jj) = 0._wp ! left side 170 IF( at_i(ji,jj+1) == 0._wp ) v_ice(ji,jj ) = 0._wp ! upper side 171 IF( at_i(ji,jj-1) == 0._wp ) v_ice(ji,jj-1) = 0._wp ! bottom side 172 ENDIF 173 END_2D 174 CALL lbc_lnk( 'icethd', u_ice, 'U', -1.0_wp, v_ice, 'V', -1.0_wp ) 279 175 ! 280 176 ! convergence tests … … 355 251 END SUBROUTINE ice_thd_mono 356 252 357 358 253 SUBROUTINE ice_thd_1d2d( kl, kn ) 359 254 !!----------------------------------------------------------------------- … … 536 431 CALL tab_1d_2d( npti, nptidx(1:npti), qcn_ice_bot_1d(1:npti), qcn_ice_bot(:,:,kl) ) 537 432 CALL tab_1d_2d( npti, nptidx(1:npti), qcn_ice_top_1d(1:npti), qcn_ice_top(:,:,kl) ) 433 CALL tab_1d_2d( npti, nptidx(1:npti), qml_ice_1d (1:npti), qml_ice (:,:,kl) ) 538 434 ! extensive variables 539 435 CALL tab_1d_2d( npti, nptidx(1:npti), v_i_1d (1:npti), v_i (:,:,kl) )
Note: See TracChangeset
for help on using the changeset viewer.
