Changeset 4990 for trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
- Timestamp:
- 2014-12-15T17:42:49+01:00 (9 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90 (modified) (23 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
r4872 r4990 22 22 USE phycst ! physical constants 23 23 USE dom_oce ! ocean space and time domain variables 24 USE oce , ONLY : iatte, oatte24 USE oce , ONLY : fraqsr_1lev 25 25 USE ice ! LIM: sea-ice variables 26 26 USE par_ice ! LIM: sea-ice parameters … … 43 43 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 44 44 USE timing ! Timing 45 USE cpl_oasis3, ONLY : lk_cpl46 45 USE limcons ! conservation tests 47 46 … … 51 50 PUBLIC lim_thd ! called by limstp module 52 51 PUBLIC lim_thd_init ! called by iceini module 53 54 REAL(wp) :: epsi10 = 1.e-10_wp !55 52 56 53 !! * Substitutions … … 68 65 !! *** ROUTINE lim_thd *** 69 66 !! 70 !! ** Purpose : This routine manages the ice thermodynamic.67 !! ** Purpose : This routine manages ice thermodynamics 71 68 !! 72 69 !! ** Action : - Initialisation of some variables … … 74 71 !! at the ice base, snow acc.,heat budget of the leads) 75 72 !! - selection of the icy points and put them in an array 76 !! - call lim_vert_ther for vert ice thermodynamic 77 !! - back to the geographic grid 78 !! - selection of points for lateral accretion 79 !! - call lim_lat_acc for the ice accretion 73 !! - call lim_thd_dif for vertical heat diffusion 74 !! - call lim_thd_dh for vertical ice growth and melt 75 !! - call lim_thd_ent for enthalpy remapping 76 !! - call lim_thd_sal for ice desalination 77 !! - call lim_thd_temp to retrieve temperature from ice enthalpy 80 78 !! - back to the geographic grid 81 79 !! 82 !! ** References : H. Goosse et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-9080 !! ** References : 83 81 !!--------------------------------------------------------------------- 84 82 INTEGER, INTENT(in) :: kt ! number of iteration … … 89 87 REAL(wp) :: zfric_umin = 0._wp ! lower bound for the friction velocity (cice value=5.e-04) 90 88 REAL(wp) :: zch = 0.0057_wp ! heat transfer coefficient 91 REAL(wp) :: z inda, zindb, zareamin89 REAL(wp) :: zareamin 92 90 REAL(wp) :: zfric_u, zqld, zqfr 93 91 ! 94 92 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 93 ! 94 REAL(wp), POINTER, DIMENSION(:,:) :: zqsr, zqns 95 95 !!------------------------------------------------------------------- 96 CALL wrk_alloc( jpi, jpj, zqsr, zqns ) 97 96 98 IF( nn_timing == 1 ) CALL timing_start('limthd') 97 99 … … 99 101 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 100 102 101 !------------------------------------------------------------------------------! 102 ! 1) Initialization of diagnostic variables ! 103 !------------------------------------------------------------------------------! 103 !------------------------------------------------------------------------! 104 ! 1) Initialization of some variables ! 105 !------------------------------------------------------------------------! 106 ftr_ice(:,:,:) = 0._wp ! part of solar radiation transmitted through the ice 107 104 108 105 109 !-------------------- … … 112 116 DO ji = 1, jpi 113 117 !0 if no ice and 1 if yes 114 zindb= 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi10 ) )118 rswitch = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi10 ) ) 115 119 !Energy of melting q(S,T) [J.m-3] 116 e_i(ji,jj,jk,jl) = zindb* e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi10 ) ) * REAL( nlay_i )120 e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi10 ) ) * REAL( nlay_i ) 117 121 !convert units ! very important that this line is here 118 122 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac … … 124 128 DO ji = 1, jpi 125 129 !0 if no ice and 1 if yes 126 zindb= 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_s(ji,jj,jl) + epsi10 ) )130 rswitch = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_s(ji,jj,jl) + epsi10 ) ) 127 131 !Energy of melting q(S,T) [J.m-3] 128 e_s(ji,jj,jk,jl) = zindb* e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi10 ) ) * REAL( nlay_s )132 e_s(ji,jj,jk,jl) = rswitch * e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi10 ) ) * REAL( nlay_s ) 129 133 !convert units ! very important that this line is here 130 134 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * unit_fac … … 136 140 ! 2) Partial computation of forcing for the thermodynamic sea ice model. ! 137 141 !-----------------------------------------------------------------------------! 142 143 !--- Ocean solar and non solar fluxes to be used in zqld 144 IF ( .NOT. lk_cpl ) THEN ! --- forced case, fluxes to the lead are the same as over the ocean 145 ! 146 zqsr(:,:) = qsr(:,:) ; zqns(:,:) = qns(:,:) 147 ! 148 ELSE ! --- coupled case, fluxes to the lead are total - intercepted 149 ! 150 zqsr(:,:) = qsr_tot(:,:) ; zqns(:,:) = qns_tot(:,:) 151 ! 152 DO jl = 1, jpl 153 DO jj = 1, jpj 154 DO ji = 1, jpi 155 zqsr(ji,jj) = zqsr(ji,jj) - qsr_ice(ji,jj,jl) * a_i_b(ji,jj,jl) 156 zqns(ji,jj) = zqns(ji,jj) - qns_ice(ji,jj,jl) * a_i_b(ji,jj,jl) 157 END DO 158 END DO 159 END DO 160 ! 161 ENDIF 138 162 139 163 !CDIR NOVERRCHK … … 141 165 !CDIR NOVERRCHK 142 166 DO ji = 1, jpi 143 zinda= tms(ji,jj) * ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - at_i(ji,jj) + epsi10 ) ) ) ! 0 if no ice167 rswitch = tms(ji,jj) * ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - at_i(ji,jj) + epsi10 ) ) ) ! 0 if no ice 144 168 ! 145 169 ! ! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget … … 149 173 ! ! temperature and turbulent mixing (McPhee, 1992) 150 174 ! 175 151 176 ! --- Energy received in the lead, zqld is defined everywhere (J.m-2) --- ! 152 zqld = tms(ji,jj) * rdt_ice * & 153 & ( pfrld(ji,jj) * ( qsr(ji,jj) * oatte(ji,jj) & ! solar heat + clem modif 154 & + qns(ji,jj) ) & ! non solar heat 155 ! latent heat of precip (note that precip is included in qns but not in qns_ice) 156 & + ( pfrld(ji,jj)**betas - pfrld(ji,jj) ) * sprecip(ji,jj) & 157 & * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rtt ) - lfus ) & 158 & + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) & 159 & * rcp * ( tatm_ice(ji,jj) - rtt ) ) 177 ! REMARK valid at least in forced mode from clem 178 ! precip is included in qns but not in qns_ice 179 IF ( lk_cpl ) THEN 180 zqld = tms(ji,jj) * rdt_ice * & 181 & ( zqsr(ji,jj) * fraqsr_1lev(ji,jj) + zqns(ji,jj) & ! pfrld already included in coupled mode 182 & + ( pfrld(ji,jj)**betas - pfrld(ji,jj) ) * sprecip(ji,jj) * & ! heat content of precip 183 & ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rtt ) - lfus ) & 184 & + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rtt ) ) 185 ELSE 186 zqld = tms(ji,jj) * rdt_ice * & 187 & ( pfrld(ji,jj) * ( zqsr(ji,jj) * fraqsr_1lev(ji,jj) + zqns(ji,jj) ) & 188 & + ( pfrld(ji,jj)**betas - pfrld(ji,jj) ) * sprecip(ji,jj) * & ! heat content of precip 189 & ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rtt ) - lfus ) & 190 & + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rtt ) ) 191 ENDIF 160 192 161 193 !-- Energy needed to bring ocean surface layer until its freezing (<0, J.m-2) --- ! … … 169 201 fhld (ji,jj) = zqld * r1_rdtice / at_i(ji,jj) ! divided by at_i since this is (re)multiplied by a_i in limthd_dh.F90 170 202 qlead(ji,jj) = 0._wp 203 ELSE 204 fhld (ji,jj) = 0._wp 171 205 ENDIF 172 206 ! … … 174 208 !clem zfric_u = MAX ( MIN( SQRT( ust2s(ji,jj) ) , zfric_umax ) , zfric_umin ) 175 209 zfric_u = MAX( SQRT( ust2s(ji,jj) ), zfric_umin ) 176 fhtur(ji,jj) = MAX( 0._wp, zinda* rau0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ) ! W.m-2210 fhtur(ji,jj) = MAX( 0._wp, rswitch * rau0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ) ! W.m-2 177 211 ! upper bound for fhtur: we do not want SST to drop below Tfreeze. 178 212 ! So we say that the heat retrieved from the ocean (fhtur+fhld) must be < to the heat necessary to reach Tfreeze (zqfr) 179 213 ! This is not a clean budget, so that should be corrected at some point 180 fhtur(ji,jj) = zinda* MIN( fhtur(ji,jj), - fhld(ji,jj) - zqfr * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) )214 fhtur(ji,jj) = rswitch * MIN( fhtur(ji,jj), - fhld(ji,jj) - zqfr * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) ) 181 215 182 216 ! ----------------------------------------- … … 187 221 hfx_in(ji,jj) = hfx_in(ji,jj) & 188 222 ! heat flux above the ocean 189 & + pfrld(ji,jj) * ( qns(ji,jj) + qsr(ji,jj) )&223 & + pfrld(ji,jj) * ( zqns(ji,jj) + zqsr(ji,jj) ) & 190 224 ! latent heat of precip (note that precip is included in qns but not in qns_ice) 191 225 & + ( 1._wp - pfrld(ji,jj) ) * sprecip(ji,jj) * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rtt ) - lfus ) & … … 198 232 ! Second step in limthd_dh : heat remaining if total melt (zq_rema) 199 233 ! Third step in limsbc : heat from ice-ocean mass exchange (zf_mass) + solar 200 hfx_out(ji,jj) = hfx_out(ji,jj) &234 hfx_out(ji,jj) = hfx_out(ji,jj) & 201 235 ! Non solar heat flux received by the ocean 202 & + pfrld(ji,jj) * qns(ji,jj) &236 & + pfrld(ji,jj) * qns(ji,jj) & 203 237 ! latent heat of precip (note that precip is included in qns but not in qns_ice) 204 & + ( pfrld(ji,jj)**betas - pfrld(ji,jj) ) * sprecip(ji,jj) 205 & * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rtt ) - lfus )&206 & + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rtt ) &238 & + ( pfrld(ji,jj)**betas - pfrld(ji,jj) ) * sprecip(ji,jj) & 239 & * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rtt ) - lfus ) & 240 & + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rtt ) & 207 241 ! heat flux taken from the ocean where there is open water ice formation 208 & - qlead(ji,jj) * r1_rdtice &242 & - qlead(ji,jj) * r1_rdtice & 209 243 ! heat flux taken from the ocean during bottom growth/melt (fhld should be 0 while bott growth) 210 & - at_i(ji,jj) * fhtur(ji,jj) &244 & - at_i(ji,jj) * fhtur(ji,jj) & 211 245 & - at_i(ji,jj) * fhld(ji,jj) 212 246 … … 310 344 CALL tab_2d_1d( nbpb, sfx_res_1d (1:nbpb), sfx_res , jpi, jpj, npb(1:nbpb) ) 311 345 312 CALL tab_2d_1d( nbpb, iatte_1d (1:nbpb), iatte , jpi, jpj, npb(1:nbpb) )313 CALL tab_2d_1d( nbpb, oatte_1d (1:nbpb), oatte , jpi, jpj, npb(1:nbpb) )314 315 346 CALL tab_2d_1d( nbpb, hfx_thd_1d (1:nbpb), hfx_thd , jpi, jpj, npb(1:nbpb) ) 316 347 CALL tab_2d_1d( nbpb, hfx_spr_1d (1:nbpb), hfx_spr , jpi, jpj, npb(1:nbpb) ) … … 389 420 CALL tab_1d_2d( nbpb, sfx_sni , npb, sfx_sni_1d(1:nbpb) , jpi, jpj ) 390 421 CALL tab_1d_2d( nbpb, sfx_res , npb, sfx_res_1d(1:nbpb) , jpi, jpj ) 391 !392 IF( num_sal == 2 ) THEN393 422 CALL tab_1d_2d( nbpb, sfx_bri , npb, sfx_bri_1d(1:nbpb) , jpi, jpj ) 394 ENDIF395 423 396 424 CALL tab_1d_2d( nbpb, hfx_thd , npb, hfx_thd_1d(1:nbpb) , jpi, jpj ) … … 407 435 CALL tab_1d_2d( nbpb, hfx_err_rem , npb, hfx_err_rem_1d(1:nbpb) , jpi, jpj ) 408 436 ! 409 !+++++ temporary stuff for a dummy version410 CALL tab_1d_2d( nbpb, dh_i_surf2D, npb, dh_i_surf(1:nbpb) , jpi, jpj )411 CALL tab_1d_2d( nbpb, dh_i_bott2D, npb, dh_i_bott(1:nbpb) , jpi, jpj )412 CALL tab_1d_2d( nbpb, s_i_newice , npb, s_i_new (1:nbpb) , jpi, jpj )413 CALL tab_1d_2d( nbpb, izero(:,:,jl) , npb, i0 (1:nbpb) , jpi, jpj )414 !+++++415 437 CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qns_ice_1d(1:nbpb) , jpi, jpj) 416 438 CALL tab_1d_2d( nbpb, ftr_ice(:,:,jl), npb, ftr_ice_1d(1:nbpb) , jpi, jpj ) … … 485 507 ENDIF 486 508 ! 509 ! 510 CALL wrk_dealloc( jpi, jpj, zqsr, zqns ) 511 512 ! 487 513 ! conservation test 488 514 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 489 515 ! 490 516 IF( nn_timing == 1 ) CALL timing_stop('limthd') 517 491 518 END SUBROUTINE lim_thd 492 519 … … 502 529 !! 503 530 INTEGER :: ji, jk ! dummy loop indices 504 REAL(wp) :: ztmelts, z switch, zaaa, zbbb, zccc, zdiscrim ! local scalar531 REAL(wp) :: ztmelts, zaaa, zbbb, zccc, zdiscrim ! local scalar 505 532 !!------------------------------------------------------------------- 506 533 ! Recover ice temperature … … 513 540 zccc = lfus * ( ztmelts - rtt ) 514 541 zdiscrim = SQRT( MAX( zbbb * zbbb - 4._wp * zaaa * zccc, 0._wp ) ) 515 t_i_1d(ji,jk) 542 t_i_1d(ji,jk) = rtt - ( zbbb + zdiscrim ) / ( 2._wp * zaaa ) 516 543 517 544 ! mask temperature 518 zswitch= 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) )519 t_i_1d(ji,jk) = zswitch * t_i_1d(ji,jk) + ( 1._wp - zswitch ) * rtt545 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) ) 546 t_i_1d(ji,jk) = rswitch * t_i_1d(ji,jk) + ( 1._wp - rswitch ) * rtt 520 547 END DO 521 548 END DO … … 555 582 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicethd in configuration namelist', lwp ) 556 583 IF(lwm) WRITE ( numoni, namicethd ) 584 585 IF( lk_cpl .AND. parsub /= 0.0 ) CALL ctl_stop( 'In coupled mode, use parsub = 0. or send dqla' ) 557 586 ! 558 587 IF(lwp) THEN ! control print
Note: See TracChangeset
for help on using the changeset viewer.
