- Timestamp:
- 2013-12-11T18:34:00+01:00 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
r4295 r4333 50 50 PUBLIC lim_thd_init ! called by iceini module 51 51 52 REAL(wp) :: epsi20 = 1e-20_wp ! constant values 53 REAL(wp) :: epsi16 = 1e-16_wp ! 54 REAL(wp) :: epsi10 = 1e-10_wp ! 55 REAL(wp) :: epsi06 = 1e-06_wp ! 56 REAL(wp) :: epsi04 = 1e-04_wp ! 52 REAL(wp) :: epsi10 = 1.e-10_wp ! 57 53 REAL(wp) :: zzero = 0._wp ! 58 54 REAL(wp) :: zone = 1._wp ! … … 126 122 DO ji = 1, jpi 127 123 !Energy of melting q(S,T) [J.m-3] 128 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi 06) ) * REAL( nlay_i )124 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi10 ) ) * REAL( nlay_i ) 129 125 !0 if no ice and 1 if yes 130 zindb = 1.0 - MAX( 0.0 , SIGN( 1.0 , - ht_i(ji,jj,jl) ) )126 zindb = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi10 ) ) 131 127 !convert units ! very important that this line is here 132 128 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac * zindb … … 138 134 DO ji = 1, jpi 139 135 !Energy of melting q(S,T) [J.m-3] 140 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi 06) ) * REAL( nlay_s )136 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi10 ) ) * REAL( nlay_s ) 141 137 !0 if no ice and 1 if yes 142 zindb = 1.0 - MAX( 0.0 , SIGN( 1.0 , - ht_s(ji,jj,jl) ) )138 zindb = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_s(ji,jj,jl) + epsi10 ) ) 143 139 !convert units ! very important that this line is here 144 140 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * unit_fac * zindb … … 147 143 END DO 148 144 END DO 149 150 !-----------------------------151 ! 1.3) Set some dummies to 0152 !-----------------------------153 !clem rdvosif(:,:) = 0.e0 ! variation of ice volume at surface154 !clem rdvobif(:,:) = 0.e0 ! variation of ice volume at bottom155 !clem fdvolif(:,:) = 0.e0 ! total variation of ice volume156 !clem rdvonif(:,:) = 0.e0 ! lateral variation of ice volume157 !clem fstric (:,:) = 0.e0 ! part of solar radiation transmitted through the ice158 !clem ffltbif(:,:) = 0.e0 ! linked with fstric159 !clem qfvbq (:,:) = 0.e0 ! linked with fstric160 !clem rdm_snw(:,:) = 0.e0 ! variation of snow mass per unit area161 !clem rdm_ice(:,:) = 0.e0 ! variation of ice mass per unit area162 !clem hicifp (:,:) = 0.e0 ! daily thermodynamic ice production.163 !clem sfx_bri(:,:) = 0.e0 ! brine flux contribution to salt flux to the ocean164 !clem fhbri (:,:) = 0.e0 ! brine flux contribution to heat flux to the ocean165 !clem sfx_thd(:,:) = 0.e0 ! equivalent salt flux to the ocean due to ice/growth decay166 145 167 146 !----------------------------------- … … 182 161 !CDIR NOVERRCHK 183 162 DO ji = 1, jpi 184 phicif(ji,jj) = vt_i(ji,jj) 185 pfrld(ji,jj) = 1.0 - at_i(ji,jj) 186 zinda = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - at_i(ji,jj) ) ) ) 163 zinda = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - at_i(ji,jj) + epsi10 ) ) ) 187 164 ! 188 165 ! ! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget … … 195 172 196 173 ! here the drag will depend on ice thickness and type (0.006) 197 fdtcn(ji,jj) = zinda * rau0 * rcp * 0.006 * zfric_u * ( ( sst_m(ji,jj) + rt0) - t_bo(ji,jj) )174 fdtcn(ji,jj) = zinda * rau0 * rcp * 0.006 * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) 198 175 ! also category dependent 199 176 ! !-- Energy from the turbulent oceanic heat flux heat flux coming in the lead 200 qdtcn(ji,jj) = zinda * fdtcn(ji,jj) * ( 1.0 - at_i(ji,jj)) * rdt_ice177 qdtcn(ji,jj) = zinda * fdtcn(ji,jj) * ( 1.0 - at_i(ji,jj) ) * rdt_ice 201 178 ! 202 179 ! !-- Lead heat budget, qldif (part 1, next one is in limthd_dh) … … 214 191 zpareff = 1.0 - zinda * zfntlat 215 192 != 0 if ice and positive heat budget and 1 if one of those two is false 216 zqlbsbq(ji,jj) = qldif(ji,jj) * ( 1.0 - zpareff ) / MAX( at_i(ji,jj) * rdt_ice , epsi16)193 zqlbsbq(ji,jj) = qldif(ji,jj) * ( 1.0 - zpareff ) / ( rdt_ice * MAX( at_i(ji,jj), epsi10 ) ) 217 194 ! 218 195 ! Heat budget of the lead, energy transferred from ice to ocean … … 221 198 ! 222 199 ! Energy needed to bring ocean surface layer until its freezing (qcmif, limflx) 223 qcmif (ji,jj) = rau0 * rcp * fse3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0) )200 qcmif (ji,jj) = rau0 * rcp * fse3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) 224 201 ! 225 202 ! oceanic heat flux (limthd_dh) 226 fbif (ji,jj) = zinda * ( fsbbq(ji,jj) / MAX( at_i(ji,jj) , epsi 20 ) + fdtcn(ji,jj) )203 fbif (ji,jj) = zinda * ( fsbbq(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) + fdtcn(ji,jj) ) 227 204 ! 228 205 END DO … … 240 217 ENDIF 241 218 242 zareamin = 1.e-10219 zareamin = epsi10 243 220 nbpb = 0 244 221 DO jj = 1, jpj … … 248 225 npb(nbpb) = (jj - 1) * jpi + ji 249 226 ENDIF 250 ! debug point to follow 251 IF ( (ji.eq.jiindx).AND.(jj.eq.jjindx) ) THEN 252 jiindex_1d = nbpb 253 ENDIF 254 END DO 255 END DO 227 END DO 228 END DO 229 230 ! debug point to follow 231 jiindex_1d = 0 232 IF( ln_nicep ) THEN 233 DO ji = mi0(jiindx), mi1(jiindx) 234 DO jj = mj0(jjindx), mj1(jjindx) 235 jiindex_1d = (jj - 1) * jpi + ji 236 END DO 237 END DO 238 ENDIF 256 239 257 240 !------------------------------------------------------------------------------! … … 315 298 !-------------------------------- 316 299 317 IF( con_i ) CALL lim_thd_enmelt( 1, nbpb ) ! computes sea ice energy of melting318 IF( con_i ) CALL lim_thd_glohec( qt_i_in, qt_s_in, q_i_layer_in, 1, nbpb, jl )300 IF( con_i .AND. jiindex_1d > 0 ) CALL lim_thd_enmelt( 1, nbpb ) ! computes sea ice energy of melting 301 IF( con_i .AND. jiindex_1d > 0 ) CALL lim_thd_glohec( qt_i_in, qt_s_in, q_i_layer_in, 1, nbpb, jl ) 319 302 320 303 ! !---------------------------------! … … 324 307 CALL lim_thd_enmelt( 1, nbpb ) ! computes sea ice energy of melting compulsory for limthd_dh 325 308 326 IF( con_i ) CALL lim_thd_glohec ( qt_i_fin, qt_s_fin, q_i_layer_fin, 1, nbpb, jl )327 IF( con_i ) CALL lim_thd_con_dif( 1 , nbpb , jl )309 IF( con_i .AND. jiindex_1d > 0 ) CALL lim_thd_glohec ( qt_i_fin, qt_s_fin, q_i_layer_fin, 1, nbpb, jl ) 310 IF( con_i .AND. jiindex_1d > 0 ) CALL lim_thd_con_dif( 1 , nbpb , jl ) 328 311 329 312 ! !---------------------------------! … … 340 323 341 324 ! CALL lim_thd_enmelt(1,nbpb) ! computes sea ice energy of melting 342 IF( con_i ) CALL lim_thd_glohec( qt_i_fin, qt_s_fin, q_i_layer_fin, 1, nbpb, jl )343 IF( con_i ) CALL lim_thd_con_dh ( 1 , nbpb , jl )325 IF( con_i .AND. jiindex_1d > 0 ) CALL lim_thd_glohec( qt_i_fin, qt_s_fin, q_i_layer_fin, 1, nbpb, jl ) 326 IF( con_i .AND. jiindex_1d > 0 ) CALL lim_thd_con_dh ( 1 , nbpb , jl ) 344 327 345 328 !-------------------------------- … … 431 414 !clem@mv-to-itd dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday 432 415 433 IF( con_i ) fbif(:,:) = fbif(:,:) + zqlbsbq(:,:)416 IF( con_i .AND. jiindex_1d > 0 ) fbif(:,:) = fbif(:,:) + zqlbsbq(:,:) 434 417 435 418 IF(ln_ctl) THEN ! Control print … … 522 505 END DO 523 506 ! 524 IF(lwp)WRITE(numout,*) ' lim_thd_glohec '525 IF(lwp)WRITE(numout,*) ' qt_i_in : ', eti(jiindex_1d,jl) * r1_rdtice526 IF(lwp)WRITE(numout,*) ' qt_s_in : ', ets(jiindex_1d,jl) * r1_rdtice527 IF(lwp)WRITE(numout,*) ' qt_in : ', ( eti(jiindex_1d,jl) + ets(jiindex_1d,jl) ) * r1_rdtice507 WRITE(numout,*) ' lim_thd_glohec ' 508 WRITE(numout,*) ' qt_i_in : ', eti(jiindex_1d,jl) * r1_rdtice 509 WRITE(numout,*) ' qt_s_in : ', ets(jiindex_1d,jl) * r1_rdtice 510 WRITE(numout,*) ' qt_in : ', ( eti(jiindex_1d,jl) + ets(jiindex_1d,jl) ) * r1_rdtice 528 511 ! 529 512 END SUBROUTINE lim_thd_glohec … … 611 594 WRITE(numout,*) ' Number of points where there is a surf err gt than surf_err : ', numce, numit 612 595 613 IF (jiindex_1D.GT.0)WRITE(numout,*) ' fc_su : ', fc_su(jiindex_1d)614 IF (jiindex_1D.GT.0)WRITE(numout,*) ' fatm : ', fatm(jiindex_1d,jl)615 IF (jiindex_1D.GT.0)WRITE(numout,*) ' t_su : ', t_su_b(jiindex_1d)596 WRITE(numout,*) ' fc_su : ', fc_su(jiindex_1d) 597 WRITE(numout,*) ' fatm : ', fatm(jiindex_1d,jl) 598 WRITE(numout,*) ' t_su : ', t_su_b(jiindex_1d) 616 599 617 600 !---------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
