Changeset 5208 for branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
- Timestamp:
- 2015-04-13T15:08:59+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
r4688 r5208 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) * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rtt ) - lfus ) & 157 & + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * 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 158 192 159 193 !-- Energy needed to bring ocean surface layer until its freezing (<0, J.m-2) --- ! … … 167 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 168 202 qlead(ji,jj) = 0._wp 203 ELSE 204 fhld (ji,jj) = 0._wp 169 205 ENDIF 170 206 ! … … 172 208 !clem zfric_u = MAX ( MIN( SQRT( ust2s(ji,jj) ) , zfric_umax ) , zfric_umin ) 173 209 zfric_u = MAX( SQRT( ust2s(ji,jj) ), zfric_umin ) 174 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 175 211 ! upper bound for fhtur: we do not want SST to drop below Tfreeze. 176 212 ! So we say that the heat retrieved from the ocean (fhtur+fhld) must be < to the heat necessary to reach Tfreeze (zqfr) 177 213 ! This is not a clean budget, so that should be corrected at some point 178 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 ) ) 179 215 180 216 ! ----------------------------------------- … … 185 221 hfx_in(ji,jj) = hfx_in(ji,jj) & 186 222 ! heat flux above the ocean 187 & + pfrld(ji,jj) * ( qns(ji,jj) + qsr(ji,jj) )&223 & + pfrld(ji,jj) * ( zqns(ji,jj) + zqsr(ji,jj) ) & 188 224 ! latent heat of precip (note that precip is included in qns but not in qns_ice) 189 225 & + ( 1._wp - pfrld(ji,jj) ) * sprecip(ji,jj) * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rtt ) - lfus ) & … … 196 232 ! Second step in limthd_dh : heat remaining if total melt (zq_rema) 197 233 ! Third step in limsbc : heat from ice-ocean mass exchange (zf_mass) + solar 198 hfx_out(ji,jj) = hfx_out(ji,jj) 234 hfx_out(ji,jj) = hfx_out(ji,jj) & 199 235 ! Non solar heat flux received by the ocean 200 & + pfrld(ji,jj) * qns(ji,jj) 236 & + pfrld(ji,jj) * qns(ji,jj) & 201 237 ! latent heat of precip (note that precip is included in qns but not in qns_ice) 202 & + ( pfrld(ji,jj)**betas - pfrld(ji,jj) ) * sprecip(ji,jj) * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rtt ) - lfus ) & 203 & + ( 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 ) & 204 241 ! heat flux taken from the ocean where there is open water ice formation 205 & - qlead(ji,jj) * r1_rdtice 242 & - qlead(ji,jj) * r1_rdtice & 206 243 ! heat flux taken from the ocean during bottom growth/melt (fhld should be 0 while bott growth) 207 & - at_i(ji,jj) * fhtur(ji,jj) 244 & - at_i(ji,jj) * fhtur(ji,jj) & 208 245 & - at_i(ji,jj) * fhld(ji,jj) 209 246 … … 256 293 !------------------------- 257 294 258 CALL tab_2d_1d( nbpb, at_i_ b(1:nbpb), at_i , jpi, jpj, npb(1:nbpb) )259 CALL tab_2d_1d( nbpb, a_i_ b(1:nbpb), a_i(:,:,jl) , jpi, jpj, npb(1:nbpb) )260 CALL tab_2d_1d( nbpb, ht_i_ b(1:nbpb), ht_i(:,:,jl) , jpi, jpj, npb(1:nbpb) )261 CALL tab_2d_1d( nbpb, ht_s_ b(1:nbpb), ht_s(:,:,jl) , jpi, jpj, npb(1:nbpb) )262 263 CALL tab_2d_1d( nbpb, t_su_ b(1:nbpb), t_su(:,:,jl) , jpi, jpj, npb(1:nbpb) )264 CALL tab_2d_1d( nbpb, sm_i_ b(1:nbpb), sm_i(:,:,jl) , jpi, jpj, npb(1:nbpb) )295 CALL tab_2d_1d( nbpb, at_i_1d (1:nbpb), at_i , jpi, jpj, npb(1:nbpb) ) 296 CALL tab_2d_1d( nbpb, a_i_1d (1:nbpb), a_i(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 297 CALL tab_2d_1d( nbpb, ht_i_1d (1:nbpb), ht_i(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 298 CALL tab_2d_1d( nbpb, ht_s_1d (1:nbpb), ht_s(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 299 300 CALL tab_2d_1d( nbpb, t_su_1d (1:nbpb), t_su(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 301 CALL tab_2d_1d( nbpb, sm_i_1d (1:nbpb), sm_i(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 265 302 DO jk = 1, nlay_s 266 CALL tab_2d_1d( nbpb, t_s_ b(1:nbpb,jk), t_s(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )267 CALL tab_2d_1d( nbpb, q_s_ b(1:nbpb,jk), e_s(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )303 CALL tab_2d_1d( nbpb, t_s_1d(1:nbpb,jk), t_s(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 304 CALL tab_2d_1d( nbpb, q_s_1d(1:nbpb,jk), e_s(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 268 305 END DO 269 306 DO jk = 1, nlay_i 270 CALL tab_2d_1d( nbpb, t_i_ b(1:nbpb,jk), t_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )271 CALL tab_2d_1d( nbpb, q_i_ b(1:nbpb,jk), e_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )272 CALL tab_2d_1d( nbpb, s_i_ b(1:nbpb,jk), s_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )307 CALL tab_2d_1d( nbpb, t_i_1d(1:nbpb,jk), t_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 308 CALL tab_2d_1d( nbpb, q_i_1d(1:nbpb,jk), e_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 309 CALL tab_2d_1d( nbpb, s_i_1d(1:nbpb,jk), s_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 273 310 END DO 274 311 … … 284 321 ENDIF 285 322 CALL tab_2d_1d( nbpb, dqns_ice_1d(1:nbpb), dqns_ice(:,:,jl), jpi, jpj, npb(1:nbpb) ) 286 CALL tab_2d_1d( nbpb, t_bo_ b(1:nbpb), t_bo , jpi, jpj, npb(1:nbpb) )323 CALL tab_2d_1d( nbpb, t_bo_1d (1:nbpb), t_bo , jpi, jpj, npb(1:nbpb) ) 287 324 CALL tab_2d_1d( nbpb, sprecip_1d (1:nbpb), sprecip , jpi, jpj, npb(1:nbpb) ) 288 325 CALL tab_2d_1d( nbpb, fhtur_1d (1:nbpb), fhtur , jpi, jpj, npb(1:nbpb) ) … … 306 343 CALL tab_2d_1d( nbpb, sfx_bri_1d (1:nbpb), sfx_bri , jpi, jpj, npb(1:nbpb) ) 307 344 CALL tab_2d_1d( nbpb, sfx_res_1d (1:nbpb), sfx_res , jpi, jpj, npb(1:nbpb) ) 308 309 CALL tab_2d_1d( nbpb, iatte_1d (1:nbpb), iatte , jpi, jpj, npb(1:nbpb) )310 CALL tab_2d_1d( nbpb, oatte_1d (1:nbpb), oatte , jpi, jpj, npb(1:nbpb) )311 345 312 346 CALL tab_2d_1d( nbpb, hfx_thd_1d (1:nbpb), hfx_thd , jpi, jpj, npb(1:nbpb) ) … … 338 372 339 373 ! --- Ice enthalpy remapping --- ! 340 CALL lim_thd_ent( 1, nbpb, q_i_ b(1:nbpb,:) )374 CALL lim_thd_ent( 1, nbpb, q_i_1d(1:nbpb,:) ) 341 375 342 376 !---------------------------------! … … 354 388 !-------------------------------- 355 389 356 CALL tab_1d_2d( nbpb, at_i , npb, at_i_ b(1:nbpb) , jpi, jpj )357 CALL tab_1d_2d( nbpb, ht_i(:,:,jl) , npb, ht_i_ b(1:nbpb) , jpi, jpj )358 CALL tab_1d_2d( nbpb, ht_s(:,:,jl) , npb, ht_s_ b(1:nbpb) , jpi, jpj )359 CALL tab_1d_2d( nbpb, a_i (:,:,jl) , npb, a_i_ b(1:nbpb) , jpi, jpj )360 CALL tab_1d_2d( nbpb, t_su(:,:,jl) , npb, t_su_ b(1:nbpb) , jpi, jpj )361 CALL tab_1d_2d( nbpb, sm_i(:,:,jl) , npb, sm_i_ b(1:nbpb) , jpi, jpj )390 CALL tab_1d_2d( nbpb, at_i , npb, at_i_1d (1:nbpb) , jpi, jpj ) 391 CALL tab_1d_2d( nbpb, ht_i(:,:,jl) , npb, ht_i_1d (1:nbpb) , jpi, jpj ) 392 CALL tab_1d_2d( nbpb, ht_s(:,:,jl) , npb, ht_s_1d (1:nbpb) , jpi, jpj ) 393 CALL tab_1d_2d( nbpb, a_i (:,:,jl) , npb, a_i_1d (1:nbpb) , jpi, jpj ) 394 CALL tab_1d_2d( nbpb, t_su(:,:,jl) , npb, t_su_1d (1:nbpb) , jpi, jpj ) 395 CALL tab_1d_2d( nbpb, sm_i(:,:,jl) , npb, sm_i_1d (1:nbpb) , jpi, jpj ) 362 396 DO jk = 1, nlay_s 363 CALL tab_1d_2d( nbpb, t_s(:,:,jk,jl), npb, t_s_ b(1:nbpb,jk), jpi, jpj)364 CALL tab_1d_2d( nbpb, e_s(:,:,jk,jl), npb, q_s_ b(1:nbpb,jk), jpi, jpj)397 CALL tab_1d_2d( nbpb, t_s(:,:,jk,jl), npb, t_s_1d (1:nbpb,jk), jpi, jpj) 398 CALL tab_1d_2d( nbpb, e_s(:,:,jk,jl), npb, q_s_1d (1:nbpb,jk), jpi, jpj) 365 399 END DO 366 400 DO jk = 1, nlay_i 367 CALL tab_1d_2d( nbpb, t_i(:,:,jk,jl), npb, t_i_ b(1:nbpb,jk), jpi, jpj)368 CALL tab_1d_2d( nbpb, e_i(:,:,jk,jl), npb, q_i_ b(1:nbpb,jk), jpi, jpj)369 CALL tab_1d_2d( nbpb, s_i(:,:,jk,jl), npb, s_i_ b(1:nbpb,jk), jpi, jpj)401 CALL tab_1d_2d( nbpb, t_i(:,:,jk,jl), npb, t_i_1d (1:nbpb,jk), jpi, jpj) 402 CALL tab_1d_2d( nbpb, e_i(:,:,jk,jl), npb, q_i_1d (1:nbpb,jk), jpi, jpj) 403 CALL tab_1d_2d( nbpb, s_i(:,:,jk,jl), npb, s_i_1d (1:nbpb,jk), jpi, jpj) 370 404 END DO 371 405 CALL tab_1d_2d( nbpb, qlead , npb, qlead_1d (1:nbpb) , jpi, jpj ) … … 386 420 CALL tab_1d_2d( nbpb, sfx_sni , npb, sfx_sni_1d(1:nbpb) , jpi, jpj ) 387 421 CALL tab_1d_2d( nbpb, sfx_res , npb, sfx_res_1d(1:nbpb) , jpi, jpj ) 388 !389 IF( num_sal == 2 ) THEN390 422 CALL tab_1d_2d( nbpb, sfx_bri , npb, sfx_bri_1d(1:nbpb) , jpi, jpj ) 391 ENDIF392 423 393 424 CALL tab_1d_2d( nbpb, hfx_thd , npb, hfx_thd_1d(1:nbpb) , jpi, jpj ) … … 404 435 CALL tab_1d_2d( nbpb, hfx_err_rem , npb, hfx_err_rem_1d(1:nbpb) , jpi, jpj ) 405 436 ! 406 !+++++ temporary stuff for a dummy version407 CALL tab_1d_2d( nbpb, dh_i_surf2D, npb, dh_i_surf(1:nbpb) , jpi, jpj )408 CALL tab_1d_2d( nbpb, dh_i_bott2D, npb, dh_i_bott(1:nbpb) , jpi, jpj )409 CALL tab_1d_2d( nbpb, s_i_newice , npb, s_i_new (1:nbpb) , jpi, jpj )410 CALL tab_1d_2d( nbpb, izero(:,:,jl) , npb, i0 (1:nbpb) , jpi, jpj )411 !+++++412 437 CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qns_ice_1d(1:nbpb) , jpi, jpj) 413 438 CALL tab_1d_2d( nbpb, ftr_ice(:,:,jl), npb, ftr_ice_1d(1:nbpb) , jpi, jpj ) … … 482 507 ENDIF 483 508 ! 509 ! 510 CALL wrk_dealloc( jpi, jpj, zqsr, zqns ) 511 512 ! 484 513 ! conservation test 485 514 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 486 515 ! 487 516 IF( nn_timing == 1 ) CALL timing_stop('limthd') 517 488 518 END SUBROUTINE lim_thd 489 519 … … 499 529 !! 500 530 INTEGER :: ji, jk ! dummy loop indices 501 REAL(wp) :: ztmelts, z switch, zaaa, zbbb, zccc, zdiscrim ! local scalar531 REAL(wp) :: ztmelts, zaaa, zbbb, zccc, zdiscrim ! local scalar 502 532 !!------------------------------------------------------------------- 503 533 ! Recover ice temperature 504 534 DO jk = 1, nlay_i 505 535 DO ji = kideb, kiut 506 ztmelts = -tmut * s_i_ b(ji,jk) + rtt536 ztmelts = -tmut * s_i_1d(ji,jk) + rtt 507 537 ! Conversion q(S,T) -> T (second order equation) 508 538 zaaa = cpic 509 zbbb = ( rcp - cpic ) * ( ztmelts - rtt ) + q_i_ b(ji,jk) / rhoic - lfus539 zbbb = ( rcp - cpic ) * ( ztmelts - rtt ) + q_i_1d(ji,jk) / rhoic - lfus 510 540 zccc = lfus * ( ztmelts - rtt ) 511 541 zdiscrim = SQRT( MAX( zbbb * zbbb - 4._wp * zaaa * zccc, 0._wp ) ) 512 t_i_ b(ji,jk)= rtt - ( zbbb + zdiscrim ) / ( 2._wp * zaaa )542 t_i_1d(ji,jk) = rtt - ( zbbb + zdiscrim ) / ( 2._wp * zaaa ) 513 543 514 544 ! mask temperature 515 zswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_i_b(ji) ) )516 t_i_ b(ji,jk) = zswitch * t_i_b(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 517 547 END DO 518 548 END DO … … 552 582 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicethd in configuration namelist', lwp ) 553 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' ) 554 586 ! 555 587 IF(lwp) THEN ! control print
Note: See TracChangeset
for help on using the changeset viewer.
