- Timestamp:
- 2014-05-12T22:46:18+02:00 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
r4332 r4634 8 8 !! 3.0 ! 2005-11 (M. Vancoppenolle) LIM-3 : Multi-layer thermodynamics + salinity variations 9 9 !! - ! 2007-04 (M. Vancoppenolle) add lim_thd_glohec, lim_thd_con_dh and lim_thd_con_dif 10 !! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in rdm_snw10 !! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in wfx_snw 11 11 !! 3.3 ! 2010-11 (G. Madec) corrected snow melting heat (due to factor betas) 12 12 !! 4.0 ! 2011-02 (G. Madec) dynamical allocation … … 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_cpl 45 46 46 47 IMPLICIT NONE … … 51 52 52 53 REAL(wp) :: epsi10 = 1.e-10_wp ! 53 REAL(wp) :: zzero = 0._wp !54 REAL(wp) :: zone = 1._wp !55 54 56 55 !! * Substitutions … … 84 83 INTEGER, INTENT(in) :: kt ! number of iteration 85 84 !! 86 INTEGER :: 87 INTEGER :: 88 REAL(wp) :: zfric_umin = 5e-03_wp ! lower bound for the friction velocity89 REAL(wp) :: zfric_umax = 2e-02_wp ! upper bound for the friction velocity90 REAL(wp) :: zinda, zindb, zthsnice, zfric_u ! local scalar91 REAL(wp) :: zfntlat, zpareff, zareamin, zcoef ! - -92 REAL(wp) , POINTER, DIMENSION(:,:) :: zqlbsbq ! link with lead energy budget qldif93 REAL(wp) :: zchk_v_i, zchk_smv, zchk_ fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset)85 INTEGER :: ji, jj, jk, jl ! dummy loop indices 86 INTEGER :: nbpb ! nb of icy pts for thermo. cal. 87 INTEGER :: ii, ij ! temporary dummy loop index 88 REAL(wp) :: zfric_umin = 5e-03_wp ! lower bound for the friction velocity 89 REAL(wp) :: zfric_umax = 2e-02_wp ! upper bound for the friction velocity 90 REAL(wp) :: zinda, zindb, zfric_u ! local scalar 91 REAL(wp) :: zareamin ! - - 92 REAL(wp) :: zchk_v_i, zchk_smv, zchk_e_i, zchk_fs, zchk_fw, zchk_ft, zchk_v_i_b, zchk_smv_b, zchk_e_i_b, zchk_fs_b, zchk_fw_b, zchk_ft_b 94 93 REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 94 REAL(wp) :: zqld, zqfr 95 REAL(wp), POINTER, DIMENSION(:) :: zdq, zq_ini, zhfx, zqfx 96 REAL(wp) :: zhfx_err, ztest 95 97 !!------------------------------------------------------------------- 96 98 IF( nn_timing == 1 ) CALL timing_start('limthd') 97 99 98 CALL wrk_alloc( jpi , jpj, zqlbsbq)100 CALL wrk_alloc( jpij, zdq, zq_ini, zhfx, zqfx ) 99 101 102 ! init debug 103 zdq(:) = 0._wp ; zq_ini(:) = 0._wp ; zhfx(:) = 0._wp ; zqfx(:) = 0._wp 104 100 105 ! ------------------------------- 101 106 !- check conservation (C Rousset) 102 107 IF (ln_limdiahsb) THEN 103 zchk_v_i_b = glob_sum( SUM( v_i(:,:,:) , dim=3 ) * area(:,:) * tms(:,:) )108 zchk_v_i_b = glob_sum( SUM( v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) 104 109 zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 105 zchk_fw_b = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) 106 zchk_fs_b = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) 110 zchk_e_i_b = glob_sum( SUM( e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) 111 zchk_fw_b = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) ) 112 zchk_fs_b = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) 113 zchk_ft_b = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) 107 114 ENDIF 108 115 !- check conservation (C Rousset) … … 121 128 DO jj = 1, jpj 122 129 DO ji = 1, jpi 123 !Energy of melting q(S,T) [J.m-3]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 )125 130 !0 if no ice and 1 if yes 126 131 zindb = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi10 ) ) 127 !convert units ! very important that this line is here 128 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac * zindb 132 !Energy of melting q(S,T) [J.m-3] 133 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 ) 134 !convert units ! very important that this line is here 135 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac 129 136 END DO 130 137 END DO … … 133 140 DO jj = 1, jpj 134 141 DO ji = 1, jpi 135 !Energy of melting q(S,T) [J.m-3]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 )137 142 !0 if no ice and 1 if yes 138 143 zindb = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_s(ji,jj,jl) + epsi10 ) ) 144 !Energy of melting q(S,T) [J.m-3] 145 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 ) 139 146 !convert units ! very important that this line is here 140 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * unit_fac * zindb147 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * unit_fac 141 148 END DO 142 149 END DO 143 150 END DO 144 151 END DO 145 146 !-----------------------------------147 ! 1.4) Compute global heat content148 !-----------------------------------149 qt_i_in (:,:) = 0.e0150 qt_s_in (:,:) = 0.e0151 qt_i_fin (:,:) = 0.e0152 qt_s_fin (:,:) = 0.e0153 sum_fluxq(:,:) = 0.e0154 fatm (:,:) = 0.e0155 152 156 153 ! 2) Partial computation of forcing for the thermodynamic sea ice model. ! … … 161 158 !CDIR NOVERRCHK 162 159 DO ji = 1, jpi 163 zinda = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - at_i(ji,jj) + epsi10 ) ) )160 zinda = tms(ji,jj) * ( 1.0 - MAX( 0._wp , SIGN( 1._wp , - at_i(ji,jj) + epsi10 ) ) ) ! 0 if no ice 164 161 ! 165 162 ! ! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget … … 168 165 ! ! net downward heat flux from the ice to the ocean, expressed as a function of ocean 169 166 ! ! temperature and turbulent mixing (McPhee, 1992) 167 170 168 ! friction velocity 171 169 zfric_u = MAX ( MIN( SQRT( ust2s(ji,jj) ) , zfric_umax ) , zfric_umin ) 172 170 173 ! here the drag will depend on ice thickness and type (0.006) 174 fdtcn(ji,jj) = zinda * rau0 * rcp * 0.006 * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) 175 ! also category dependent 176 ! !-- Energy from the turbulent oceanic heat flux heat flux coming in the lead 177 qdtcn(ji,jj) = zinda * fdtcn(ji,jj) * ( 1.0 - at_i(ji,jj) ) * rdt_ice 178 ! 179 ! !-- Lead heat budget, qldif (part 1, next one is in limthd_dh) 180 ! ! caution: exponent betas used as more snow can fallinto leads 181 qldif(ji,jj) = tms(ji,jj) * rdt_ice * ( & 182 & pfrld(ji,jj) * ( qsr(ji,jj) * oatte(ji,jj) & ! solar heat + clem modif 183 & + qns(ji,jj) & ! non solar heat 184 & + fdtcn(ji,jj) & ! turbulent ice-ocean heat 185 & + fsbbq(ji,jj) * ( 1.0 - zinda ) ) & ! residual heat from previous step 186 & - pfrld(ji,jj)**betas * sprecip(ji,jj) * lfus ) ! latent heat of sprecip melting 171 !-- Energy from the turbulent oceanic heat flux. here the drag will depend on ice thickness and type (0.006) 172 fhtur(ji,jj) = zinda * rau0 * rcp * 0.006 * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ! W.m-2 173 ! clem: why not the following? 174 !fhtur(ji,jj) = zinda * rau0 * rcp * 0.006 * SQRT( ust2s(ji,jj) ) * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) 175 176 !-- Energy received in the lead, zqld is defined everywhere (J.m-2) 177 ! It includes turbulent ocean heat flux (only in the leads, the rest is used for bottom melting) 178 zqld = tms(ji,jj) * rdt_ice * & 179 & ( pfrld(ji,jj) * ( qsr(ji,jj) * oatte(ji,jj) & ! solar heat + clem modif 180 & + qns(ji,jj) & ! non solar heat 181 & + fhtur(ji,jj) ) & ! turbulent ice-ocean heat (0 if no ice) 182 ! latent heat of precip (note that precip is included in qns but not in qns_ice) 183 & + ( pfrld(ji,jj)**betas - pfrld(ji,jj) ) * sprecip(ji,jj) * ( 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 186 !-- Energy needed to bring ocean surface layer until its freezing (<0, J.m-2) 187 zqfr = tms(ji,jj) * rau0 * rcp * fse3t_m(ji,jj,1) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) 188 189 !-- Energy Budget of the leads (J.m-2). Must be < 0 to form ice 190 qlead(ji,jj) = MIN( 0._wp , zqld - zqfr ) 191 192 ! If there is ice and leads are warming, then transfer energy from the lead budget and use it for bottom melting 193 IF( at_i(ji,jj) > epsi10 .AND. zqld > 0._wp ) THEN 194 fhld (ji,jj) = zqld * r1_rdtice / at_i(ji,jj) ! divided by a_i since this is (re)multiplied by a_i in limthd_dh.F90 195 qlead(ji,jj) = 0._wp 196 ENDIF 187 197 ! 188 ! Positive heat budget is used for bottom ablation 189 zfntlat = 1.0 - MAX( zzero , SIGN( zone , - qldif(ji,jj) ) ) 190 != 1 if positive heat budget 191 zpareff = 1.0 - zinda * zfntlat 192 != 0 if ice and positive heat budget and 1 if one of those two is false 193 zqlbsbq(ji,jj) = qldif(ji,jj) * ( 1.0 - zpareff ) / ( rdt_ice * MAX( at_i(ji,jj), epsi10 ) ) 198 IF( qlead(ji,jj) == 0._wp ) zqld = 0._wp ; zqfr = 0._wp 194 199 ! 195 ! Heat budget of the lead, energy transferred from ice to ocean 196 qldif (ji,jj) = zpareff * qldif(ji,jj) 197 qdtcn (ji,jj) = zpareff * qdtcn(ji,jj) 198 ! 199 ! Energy needed to bring ocean surface layer until its freezing (qcmif, limflx) 200 qcmif (ji,jj) = rau0 * rcp * fse3t_m(ji,jj,1) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) 201 ! 202 ! oceanic heat flux (limthd_dh) 203 fbif (ji,jj) = zinda * ( fsbbq(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) + fdtcn(ji,jj) ) 204 ! 200 ! ----------------------------------------- 201 ! Net heat flux on top of ice-ocean [W.m-2] 202 ! ----------------------------------------- 203 ! First step here : heat flux at the ocean surface + precip 204 ! Second step below : heat flux at the ice surface (after limthd_dif) 205 hfx_in(ji,jj) = hfx_in(ji,jj) & 206 ! heat flux above the ocean 207 & + pfrld(ji,jj) * ( qns(ji,jj) + qsr(ji,jj) ) & 208 ! latent heat of precip (note that precip is included in qns but not in qns_ice) 209 & + ( 1._wp - pfrld(ji,jj) ) * sprecip(ji,jj) * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rtt ) - lfus ) & 210 & + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rtt ) 211 212 ! ----------------------------------------------------------------------------- 213 ! Net heat flux that is retroceded to the ocean or taken from the ocean [W.m-2] 214 ! ----------------------------------------------------------------------------- 215 ! First step here : non solar + precip - qlead - qturb 216 ! Second step in limthd_dh : heat remaining if total melt (zq_rema) 217 ! Third step in limsbc : heat from ice-ocean mass exchange (zf_mass) + solar 218 hfx_out(ji,jj) = hfx_out(ji,jj) & 219 ! Non solar heat flux received by the ocean 220 & + pfrld(ji,jj) * qns(ji,jj) & 221 ! latent heat of precip (note that precip is included in qns but not in qns_ice) 222 & + ( pfrld(ji,jj)**betas - pfrld(ji,jj) ) * sprecip(ji,jj) * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rtt ) - lfus ) & 223 & + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rtt ) & 224 ! heat flux taken from the ocean where there is open water ice formation 225 & - qlead(ji,jj) * r1_rdtice & 226 ! heat flux taken from the ocean during bottom growth/melt (fhld should be 0 while bott growth) 227 & - at_i(ji,jj) * fhtur(ji,jj) & 228 & - at_i(ji,jj) * fhld(ji,jj) 229 205 230 END DO 206 231 END DO … … 234 259 DO jj = mj0(jjindx), mj1(jjindx) 235 260 jiindex_1d = (jj - 1) * jpi + ji 261 WRITE(numout,*) ' lim_thd : Category no : ', jl 236 262 END DO 237 263 END DO … … 271 297 CALL tab_2d_1d( nbpb, fr1_i0_1d (1:nbpb), fr1_i0 , jpi, jpj, npb(1:nbpb) ) 272 298 CALL tab_2d_1d( nbpb, fr2_i0_1d (1:nbpb), fr2_i0 , jpi, jpj, npb(1:nbpb) ) 273 CALL tab_2d_1d( nbpb, qnsr_ice_1d(1:nbpb), qns_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 274 #if ! defined key_coupled 275 CALL tab_2d_1d( nbpb, qla_ice_1d (1:nbpb), qla_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 276 CALL tab_2d_1d( nbpb, dqla_ice_1d(1:nbpb), dqla_ice(:,:,jl), jpi, jpj, npb(1:nbpb) ) 277 #endif 299 CALL tab_2d_1d( nbpb, qns_ice_1d (1:nbpb), qns_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 300 CALL tab_2d_1d( nbpb, ftr_ice_1d (1:nbpb), ftr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 301 IF( .NOT. lk_cpl ) THEN 302 CALL tab_2d_1d( nbpb, qla_ice_1d (1:nbpb), qla_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 303 CALL tab_2d_1d( nbpb, dqla_ice_1d(1:nbpb), dqla_ice(:,:,jl), jpi, jpj, npb(1:nbpb) ) 304 ENDIF 278 305 CALL tab_2d_1d( nbpb, dqns_ice_1d(1:nbpb), dqns_ice(:,:,jl), jpi, jpj, npb(1:nbpb) ) 279 306 CALL tab_2d_1d( nbpb, t_bo_b (1:nbpb), t_bo , jpi, jpj, npb(1:nbpb) ) 280 307 CALL tab_2d_1d( nbpb, sprecip_1d (1:nbpb), sprecip , jpi, jpj, npb(1:nbpb) ) 281 CALL tab_2d_1d( nbpb, fbif_1d (1:nbpb), fbif , jpi, jpj, npb(1:nbpb) ) 282 CALL tab_2d_1d( nbpb, qldif_1d (1:nbpb), qldif , jpi, jpj, npb(1:nbpb) ) 283 CALL tab_2d_1d( nbpb, rdm_ice_1d (1:nbpb), rdm_ice , jpi, jpj, npb(1:nbpb) ) 284 CALL tab_2d_1d( nbpb, rdm_snw_1d (1:nbpb), rdm_snw , jpi, jpj, npb(1:nbpb) ) 285 CALL tab_2d_1d( nbpb, dmgwi_1d (1:nbpb), dmgwi , jpi, jpj, npb(1:nbpb) ) 286 CALL tab_2d_1d( nbpb, qlbbq_1d (1:nbpb), zqlbsbq , jpi, jpj, npb(1:nbpb) ) 287 288 CALL tab_2d_1d( nbpb, sfx_thd_1d (1:nbpb), sfx_thd , jpi, jpj, npb(1:nbpb) ) 308 CALL tab_2d_1d( nbpb, fhtur_1d (1:nbpb), fhtur , jpi, jpj, npb(1:nbpb) ) 309 CALL tab_2d_1d( nbpb, qlead_1d (1:nbpb), qlead , jpi, jpj, npb(1:nbpb) ) 310 CALL tab_2d_1d( nbpb, fhld_1d (1:nbpb), fhld , jpi, jpj, npb(1:nbpb) ) 311 312 CALL tab_2d_1d( nbpb, wfx_snw_1d (1:nbpb), wfx_snw , jpi, jpj, npb(1:nbpb) ) 313 CALL tab_2d_1d( nbpb, wfx_sub_1d (1:nbpb), wfx_sub , jpi, jpj, npb(1:nbpb) ) 314 315 CALL tab_2d_1d( nbpb, wfx_bog_1d (1:nbpb), wfx_bog , jpi, jpj, npb(1:nbpb) ) 316 CALL tab_2d_1d( nbpb, wfx_bom_1d (1:nbpb), wfx_bom , jpi, jpj, npb(1:nbpb) ) 317 CALL tab_2d_1d( nbpb, wfx_sum_1d (1:nbpb), wfx_sum , jpi, jpj, npb(1:nbpb) ) 318 CALL tab_2d_1d( nbpb, wfx_sni_1d (1:nbpb), wfx_sni , jpi, jpj, npb(1:nbpb) ) 319 320 CALL tab_2d_1d( nbpb, sfx_bog_1d (1:nbpb), sfx_bog , jpi, jpj, npb(1:nbpb) ) 321 CALL tab_2d_1d( nbpb, sfx_bom_1d (1:nbpb), sfx_bom , jpi, jpj, npb(1:nbpb) ) 322 CALL tab_2d_1d( nbpb, sfx_sum_1d (1:nbpb), sfx_sum , jpi, jpj, npb(1:nbpb) ) 323 CALL tab_2d_1d( nbpb, sfx_sni_1d (1:nbpb), sfx_sni , jpi, jpj, npb(1:nbpb) ) 289 324 CALL tab_2d_1d( nbpb, sfx_bri_1d (1:nbpb), sfx_bri , jpi, jpj, npb(1:nbpb) ) 290 CALL tab_2d_1d( nbpb, fhbri_1d (1:nbpb), fhbri , jpi, jpj, npb(1:nbpb) ) 291 CALL tab_2d_1d( nbpb, fstbif_1d (1:nbpb), fstric , jpi, jpj, npb(1:nbpb) ) 292 CALL tab_2d_1d( nbpb, qfvbq_1d (1:nbpb), qfvbq , jpi, jpj, npb(1:nbpb) ) 293 294 CALL tab_2d_1d( nbpb, iatte_1d (1:nbpb), iatte , jpi, jpj, npb(1:nbpb) ) ! clem modif 295 CALL tab_2d_1d( nbpb, oatte_1d (1:nbpb), oatte , jpi, jpj, npb(1:nbpb) ) ! clem modif 325 326 CALL tab_2d_1d( nbpb, iatte_1d (1:nbpb), iatte , jpi, jpj, npb(1:nbpb) ) 327 CALL tab_2d_1d( nbpb, oatte_1d (1:nbpb), oatte , jpi, jpj, npb(1:nbpb) ) 328 329 CALL tab_2d_1d( nbpb, hfx_thd_1d (1:nbpb), hfx_thd , jpi, jpj, npb(1:nbpb) ) 330 CALL tab_2d_1d( nbpb, hfx_spr_1d (1:nbpb), hfx_spr , jpi, jpj, npb(1:nbpb) ) 331 CALL tab_2d_1d( nbpb, hfx_tot_1d (1:nbpb), hfx_tot , jpi, jpj, npb(1:nbpb) ) 332 CALL tab_2d_1d( nbpb, hfx_snw_1d (1:nbpb), hfx_snw , jpi, jpj, npb(1:nbpb) ) 333 CALL tab_2d_1d( nbpb, hfx_sub_1d (1:nbpb), hfx_sub , jpi, jpj, npb(1:nbpb) ) 334 CALL tab_2d_1d( nbpb, hfx_err_1d (1:nbpb), hfx_err , jpi, jpj, npb(1:nbpb) ) 335 CALL tab_2d_1d( nbpb, hfx_res_1d (1:nbpb), hfx_res , jpi, jpj, npb(1:nbpb) ) 336 CALL tab_2d_1d( nbpb, hfx_err_rem_1d (1:nbpb), hfx_err_rem , jpi, jpj, npb(1:nbpb) ) 337 296 338 !-------------------------------- 297 339 ! 4.3) Thermodynamic processes 298 340 !-------------------------------- 299 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 ) 302 303 ! !---------------------------------! 304 CALL lim_thd_dif( 1, nbpb, jl ) ! Ice/Snow Temperature profile ! 305 ! !---------------------------------! 306 307 CALL lim_thd_enmelt( 1, nbpb ) ! computes sea ice energy of melting compulsory for limthd_dh 308 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 ) 311 312 ! !---------------------------------! 313 CALL lim_thd_dh( 1, nbpb, jl ) ! Ice/Snow thickness ! 314 ! !---------------------------------! 315 316 ! !---------------------------------! 317 CALL lim_thd_ent( 1, nbpb, jl ) ! Ice/Snow enthalpy remapping ! 318 ! !---------------------------------! 319 320 ! !---------------------------------! 321 CALL lim_thd_sal( 1, nbpb ) ! Ice salinity computation ! 322 ! !---------------------------------! 341 ! --- diag error on heat diffusion - PART 1 --- ! 342 DO ji = 1, nbpb 343 zq_ini(ji) = ( SUM( q_i_b(ji,1:nlay_i) ) * ht_i_b(ji) / REAL( nlay_i ) + & 344 & SUM( q_s_b(ji,1:nlay_s) ) * ht_s_b(ji) / REAL( nlay_s ) ) 345 END DO 346 347 !---------------------------------! 348 ! Ice/Snow Temperature profile ! 349 !---------------------------------! 350 CALL lim_thd_dif( 1, nbpb, jl ) 351 352 ! --- computes sea ice energy of melting compulsory for limthd_dh --- ! 353 CALL lim_thd_enmelt( 1, nbpb ) 354 355 DO ji = 1, nbpb 356 ! --- diag error on heat diffusion - PART 2 --- ! 357 zdq(ji) = - zq_ini(ji) + ( SUM( q_i_b(ji,1:nlay_i) ) * ht_i_b(ji) / REAL( nlay_i ) + & 358 & SUM( q_s_b(ji,1:nlay_s) ) * ht_s_b(ji) / REAL( nlay_s ) ) 359 zhfx_err = ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - ftr_ice_1d(ji) - fc_bo_i(ji) + zdq(ji) * r1_rdtice ) 360 hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err * a_i_b(ji) 361 ! --- correction of qns_ice and surface conduction flux --- ! 362 qns_ice_1d(ji) = qns_ice_1d(ji) - zhfx_err 363 fc_su (ji) = fc_su (ji) - zhfx_err 364 ! --- Heat flux at the ice surface in W.m-2 --- ! 365 ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 366 hfx_in (ii,ij) = hfx_in (ii,ij) + a_i_b(ji) * ( qsr_ice_1d(ji) + qns_ice_1d(ji) ) 367 368 END DO 369 370 !---------------------------------! 371 ! Ice/Snow thicnkess ! 372 !---------------------------------! 373 ! --- diag error on heat remapping - PART 1 --- ! 374 DO ji = 1, nbpb 375 zq_ini(ji) = ( SUM( q_i_b(ji,1:nlay_i) ) * ht_i_b(ji) / REAL( nlay_i ) + & 376 & SUM( q_s_b(ji,1:nlay_s) ) * ht_s_b(ji) / REAL( nlay_s ) ) 377 END DO 378 379 CALL lim_thd_dh( 1, nbpb, jl ) 380 381 ! --- Ice/Snow enthalpy remapping --- ! 382 CALL lim_thd_ent( 1, nbpb, jl ) 383 ! 384 ! --- diag error on heat remapping - PART 2 --- ! 385 DO ji = 1, nbpb 386 zdq(ji) = - ( zq_ini(ji) + dq_i(ji) + dq_s(ji) ) & 387 & + ( SUM( q_i_b(ji,1:nlay_i) ) * ht_i_b(ji) / REAL( nlay_i ) + & 388 & SUM( q_s_b(ji,1:nlay_s) ) * ht_s_b(ji) / REAL( nlay_s ) ) 389 hfx_err_rem_1d(ji) = hfx_err_rem_1d(ji) + zdq(ji) * a_i_b(ji) * r1_rdtice 390 END DO 391 392 !---------------------------------! 393 ! Ice salinity ! 394 !---------------------------------! 395 CALL lim_thd_sal( 1, nbpb ) 323 396 324 397 ! CALL lim_thd_enmelt(1,nbpb) ! computes sea ice energy of melting 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 )327 328 398 !-------------------------------- 329 399 ! 4.4) Move 1D to 2D vectors … … 345 415 CALL tab_1d_2d( nbpb, s_i(:,:,jk,jl), npb, s_i_b (1:nbpb,jk), jpi, jpj) 346 416 END DO 347 CALL tab_1d_2d( nbpb, fstric , npb, fstbif_1d (1:nbpb) , jpi, jpj ) 348 CALL tab_1d_2d( nbpb, qldif , npb, qldif_1d (1:nbpb) , jpi, jpj ) 349 CALL tab_1d_2d( nbpb, qfvbq , npb, qfvbq_1d (1:nbpb) , jpi, jpj ) 350 CALL tab_1d_2d( nbpb, rdm_ice , npb, rdm_ice_1d(1:nbpb) , jpi, jpj ) 351 CALL tab_1d_2d( nbpb, rdm_snw , npb, rdm_snw_1d(1:nbpb) , jpi, jpj ) 352 CALL tab_1d_2d( nbpb, dmgwi , npb, dmgwi_1d (1:nbpb) , jpi, jpj ) 353 CALL tab_1d_2d( nbpb, rdvosif , npb, dvsbq_1d (1:nbpb) , jpi, jpj ) 354 CALL tab_1d_2d( nbpb, rdvobif , npb, dvbbq_1d (1:nbpb) , jpi, jpj ) 355 CALL tab_1d_2d( nbpb, fdvolif , npb, dvlbq_1d (1:nbpb) , jpi, jpj ) 356 CALL tab_1d_2d( nbpb, rdvonif , npb, dvnbq_1d (1:nbpb) , jpi, jpj ) 357 CALL tab_1d_2d( nbpb, sfx_thd , npb, sfx_thd_1d(1:nbpb) , jpi, jpj ) 417 CALL tab_1d_2d( nbpb, qlead , npb, qlead_1d (1:nbpb) , jpi, jpj ) 418 419 CALL tab_1d_2d( nbpb, wfx_snw , npb, wfx_snw_1d(1:nbpb) , jpi, jpj ) 420 CALL tab_1d_2d( nbpb, wfx_sub , npb, wfx_sub_1d(1:nbpb) , jpi, jpj ) 421 422 CALL tab_1d_2d( nbpb, wfx_bog , npb, wfx_bog_1d(1:nbpb) , jpi, jpj ) 423 CALL tab_1d_2d( nbpb, wfx_bom , npb, wfx_bom_1d(1:nbpb) , jpi, jpj ) 424 CALL tab_1d_2d( nbpb, wfx_sum , npb, wfx_sum_1d(1:nbpb) , jpi, jpj ) 425 CALL tab_1d_2d( nbpb, wfx_sni , npb, wfx_sni_1d(1:nbpb) , jpi, jpj ) 426 427 CALL tab_1d_2d( nbpb, sfx_bog , npb, sfx_bog_1d(1:nbpb) , jpi, jpj ) 428 CALL tab_1d_2d( nbpb, sfx_bom , npb, sfx_bom_1d(1:nbpb) , jpi, jpj ) 429 CALL tab_1d_2d( nbpb, sfx_sum , npb, sfx_sum_1d(1:nbpb) , jpi, jpj ) 430 CALL tab_1d_2d( nbpb, sfx_sni , npb, sfx_sni_1d(1:nbpb) , jpi, jpj ) 358 431 ! 359 432 IF( num_sal == 2 ) THEN 360 433 CALL tab_1d_2d( nbpb, sfx_bri , npb, sfx_bri_1d(1:nbpb) , jpi, jpj ) 361 CALL tab_1d_2d( nbpb, fhbri , npb, fhbri_1d (1:nbpb) , jpi, jpj )362 434 ENDIF 435 436 CALL tab_1d_2d( nbpb, hfx_thd , npb, hfx_thd_1d(1:nbpb) , jpi, jpj ) 437 CALL tab_1d_2d( nbpb, hfx_spr , npb, hfx_spr_1d(1:nbpb) , jpi, jpj ) 438 CALL tab_1d_2d( nbpb, hfx_tot , npb, hfx_tot_1d(1:nbpb) , jpi, jpj ) 439 CALL tab_1d_2d( nbpb, hfx_snw , npb, hfx_snw_1d(1:nbpb) , jpi, jpj ) 440 CALL tab_1d_2d( nbpb, hfx_sub , npb, hfx_sub_1d(1:nbpb) , jpi, jpj ) 441 CALL tab_1d_2d( nbpb, hfx_err , npb, hfx_err_1d(1:nbpb) , jpi, jpj ) 442 CALL tab_1d_2d( nbpb, hfx_res , npb, hfx_res_1d(1:nbpb) , jpi, jpj ) 443 CALL tab_1d_2d( nbpb, hfx_err_rem , npb, hfx_err_rem_1d(1:nbpb) , jpi, jpj ) 363 444 ! 364 445 !+++++ temporary stuff for a dummy version 365 CALL tab_1d_2d( nbpb, dh_i_surf2D, npb, dh_i_surf(1:nbpb) , jpi, jpj ) 366 CALL tab_1d_2d( nbpb, dh_i_bott2D, npb, dh_i_bott(1:nbpb) , jpi, jpj ) 367 CALL tab_1d_2d( nbpb, fsup2D , npb, fsup (1:nbpb) , jpi, jpj ) 368 CALL tab_1d_2d( nbpb, focea2D , npb, focea (1:nbpb) , jpi, jpj ) 369 CALL tab_1d_2d( nbpb, s_i_newice , npb, s_i_new (1:nbpb) , jpi, jpj ) 370 CALL tab_1d_2d( nbpb, izero(:,:,jl) , npb, i0 (1:nbpb) , jpi, jpj ) 371 CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qnsr_ice_1d(1:nbpb), jpi, jpj) 446 CALL tab_1d_2d( nbpb, dh_i_surf2D, npb, dh_i_surf(1:nbpb) , jpi, jpj ) 447 CALL tab_1d_2d( nbpb, dh_i_bott2D, npb, dh_i_bott(1:nbpb) , jpi, jpj ) 448 CALL tab_1d_2d( nbpb, s_i_newice , npb, s_i_new (1:nbpb) , jpi, jpj ) 449 CALL tab_1d_2d( nbpb, izero(:,:,jl) , npb, i0 (1:nbpb) , jpi, jpj ) 372 450 !+++++ 451 CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qns_ice_1d(1:nbpb) , jpi, jpj) 452 CALL tab_1d_2d( nbpb, ftr_ice(:,:,jl), npb, ftr_ice_1d(1:nbpb) , jpi, jpj ) 373 453 ! 374 454 IF( lk_mpp ) CALL mpp_comm_free( ncomm_ice ) !RB necessary ?? … … 384 464 ! 5.1) Ice heat content 385 465 !------------------------ 386 ! Enthalpies are global variables we have to readjust the units (heat content in 10^9 Joules) 387 zcoef = 1._wp / ( unit_fac * REAL( nlay_i ) ) 466 ! Enthalpies are global variables we have to readjust the units (heat content in Joules) 388 467 DO jl = 1, jpl 389 468 DO jk = 1, nlay_i 390 e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * a_i(:,:,jl) * ht_i(:,:,jl) * zcoef469 e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * a_i(:,:,jl) * ht_i(:,:,jl) / ( unit_fac * REAL( nlay_i ) ) 391 470 END DO 392 471 END DO … … 395 474 ! 5.2) Snow heat content 396 475 !------------------------ 397 ! Enthalpies are global variables we have to readjust the units (heat content in 10^9 Joules) 398 zcoef = 1._wp / ( unit_fac * REAL( nlay_s ) ) 476 ! Enthalpies are global variables we have to readjust the units (heat content in Joules) 399 477 DO jl = 1, jpl 400 478 DO jk = 1, nlay_s 401 e_s(:,:,jk,jl) = e_s(:,:,jk,jl) * area(:,:) * a_i(:,:,jl) * ht_s(:,:,jl) * zcoef479 e_s(:,:,jk,jl) = e_s(:,:,jk,jl) * area(:,:) * a_i(:,:,jl) * ht_s(:,:,jl) / ( unit_fac * REAL( nlay_s ) ) 402 480 END DO 403 481 END DO … … 411 489 ! 5.4) Diagnostic thermodynamic growth rates 412 490 !-------------------------------------------- 413 !clem@useless d_v_i_thd(:,:,:) = v_i (:,:,:) - old_v_i(:,:,:) ! ice volumes414 !clem@mv-to-itd dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday415 416 IF( con_i .AND. jiindex_1d > 0 ) fbif(:,:) = fbif(:,:) + zqlbsbq(:,:)417 418 491 IF(ln_ctl) THEN ! Control print 419 492 CALL prt_ctl_info(' ') … … 451 524 !- check conservation (C Rousset) 452 525 IF (ln_limdiahsb) THEN 453 zchk_fs = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 454 zchk_fw = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b 526 zchk_fs = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 527 zchk_fw = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fw_b 528 zchk_ft = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) - zchk_ft_b 455 529 456 zchk_v_i = ( glob_sum( SUM( v_i(:,:,:) , dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) * r1_rdtice530 zchk_v_i = ( glob_sum( SUM( v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b ) * r1_rdtice - zchk_fw 457 531 zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic ) 532 zchk_e_i = glob_sum( SUM( e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) * r1_rdtice - zchk_e_i_b * r1_rdtice + zchk_ft 458 533 459 534 zchk_vmin = glob_min(v_i) … … 462 537 463 538 IF(lwp) THEN 464 IF ( ABS( zchk_v_i ) > 1.e- 5 ) WRITE(numout,*) 'violation volume [m3/day] (limthd) = ',(zchk_v_i * rday)539 IF ( ABS( zchk_v_i ) > 1.e-4 ) WRITE(numout,*) 'violation volume [kg/day] (limthd) = ',(zchk_v_i * rday) 465 540 IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limthd) = ',(zchk_smv * rday) 541 IF ( ABS( zchk_e_i ) > 1.e-2 ) WRITE(numout,*) 'violation enthalpy [1e9 J] (limthd) = ',(zchk_e_i) 466 542 IF ( zchk_vmin < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limthd) = ',(zchk_vmin * 1.e-3) 467 543 IF ( zchk_amax > amax+epsi10 ) WRITE(numout,*) 'violation a_i>amax (limthd) = ',zchk_amax … … 472 548 ! ------------------------------- 473 549 ! 474 CALL wrk_dealloc( jpi , jpj, zqlbsbq)475 ! 550 CALL wrk_dealloc( jpij, zdq, zq_ini, zhfx, zqfx ) 551 476 552 IF( nn_timing == 1 ) CALL timing_stop('limthd') 477 553 END SUBROUTINE lim_thd 478 554 479 480 SUBROUTINE lim_thd_glohec( eti, ets, etilayer, kideb, kiut, jl ) 481 !!----------------------------------------------------------------------- 482 !! *** ROUTINE lim_thd_glohec *** 483 !! 484 !! ** Purpose : Compute total heat content for each category 485 !! Works with 1d vectors only 486 !!----------------------------------------------------------------------- 487 INTEGER , INTENT(in ) :: kideb, kiut ! bounds for the spatial loop 488 INTEGER , INTENT(in ) :: jl ! category number 489 REAL(wp), INTENT( out), DIMENSION (jpij,jpl ) :: eti, ets ! vertically-summed heat content for ice & snow 490 REAL(wp), INTENT( out), DIMENSION (jpij,jkmax) :: etilayer ! heat content for ice layers 491 !! 492 INTEGER :: ji,jk ! loop indices 493 !!----------------------------------------------------------------------- 494 eti(:,:) = 0._wp 495 ets(:,:) = 0._wp 496 ! 497 DO jk = 1, nlay_i ! total q over all layers, ice [J.m-2] 498 DO ji = kideb, kiut 499 etilayer(ji,jk) = q_i_b(ji,jk) * ht_i_b(ji) / REAL( nlay_i ) 500 eti (ji,jl) = eti(ji,jl) + etilayer(ji,jk) 501 END DO 502 END DO 503 DO ji = kideb, kiut ! total q over all layers, snow [J.m-2] 504 ets(ji,jl) = ets(ji,jl) + q_s_b(ji,1) * ht_s_b(ji) / REAL( nlay_s ) 505 END DO 506 ! 507 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 511 ! 512 END SUBROUTINE lim_thd_glohec 513 514 515 SUBROUTINE lim_thd_con_dif( kideb, kiut, jl ) 516 !!----------------------------------------------------------------------- 517 !! *** ROUTINE lim_thd_con_dif *** 518 !! 519 !! ** Purpose : Test energy conservation after heat diffusion 520 !!------------------------------------------------------------------- 521 INTEGER , INTENT(in ) :: kideb, kiut ! bounds for the spatial loop 522 INTEGER , INTENT(in ) :: jl ! category number 523 524 INTEGER :: ji, jk ! loop indices 525 INTEGER :: ii, ij 526 INTEGER :: numce ! number of points for which conservation is violated 527 REAL(wp) :: meance ! mean conservation error 528 REAL(wp) :: max_cons_err, max_surf_err 529 !!--------------------------------------------------------------------- 530 531 max_cons_err = 1.0_wp ! maximum tolerated conservation error 532 max_surf_err = 0.001_wp ! maximum tolerated surface error 533 534 !-------------------------- 535 ! Increment of energy 536 !-------------------------- 537 ! global 538 DO ji = kideb, kiut 539 dq_i(ji,jl) = qt_i_fin(ji,jl) - qt_i_in(ji,jl) + qt_s_fin(ji,jl) - qt_s_in(ji,jl) 540 END DO 541 ! layer by layer 542 dq_i_layer(:,:) = q_i_layer_fin(:,:) - q_i_layer_in(:,:) 543 544 !---------------------------------------- 545 ! Atmospheric heat flux, ice heat budget 546 !---------------------------------------- 547 DO ji = kideb, kiut 548 ii = MOD( npb(ji) - 1 , jpi ) + 1 549 ij = ( npb(ji) - 1 ) / jpi + 1 550 fatm (ji,jl) = qnsr_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) 551 sum_fluxq(ji,jl) = fc_su(ji) - fc_bo_i(ji) + qsr_ice_1d(ji) * i0(ji) - fstroc(ii,ij,jl) 552 END DO 553 554 !-------------------- 555 ! Conservation error 556 !-------------------- 557 DO ji = kideb, kiut 558 cons_error(ji,jl) = ABS( dq_i(ji,jl) * r1_rdtice + sum_fluxq(ji,jl) ) 559 END DO 560 561 numce = 0 562 meance = 0._wp 563 DO ji = kideb, kiut 564 IF ( cons_error(ji,jl) .GT. max_cons_err ) THEN 565 numce = numce + 1 566 meance = meance + cons_error(ji,jl) 567 ENDIF 568 END DO 569 IF( numce > 0 ) meance = meance / numce 570 571 WRITE(numout,*) ' Maximum tolerated conservation error : ', max_cons_err 572 WRITE(numout,*) ' After lim_thd_dif, category : ', jl 573 WRITE(numout,*) ' Mean conservation error on big error points ', meance, numit 574 WRITE(numout,*) ' Number of points where there is a cons err gt than c.e. : ', numce, numit 575 576 !------------------------------------------------------- 577 ! Surface error due to imbalance between Fatm and Fcsu 578 !------------------------------------------------------- 579 numce = 0 580 meance = 0._wp 581 582 DO ji = kideb, kiut 583 surf_error(ji,jl) = ABS ( fatm(ji,jl) - fc_su(ji) ) 584 IF( ( t_su_b(ji) .LT. rtt ) .AND. ( surf_error(ji,jl) .GT. max_surf_err ) ) THEN 585 numce = numce + 1 586 meance = meance + surf_error(ji,jl) 587 ENDIF 588 ENDDO 589 IF( numce > 0 ) meance = meance / numce 590 591 WRITE(numout,*) ' Maximum tolerated surface error : ', max_surf_err 592 WRITE(numout,*) ' After lim_thd_dif, category : ', jl 593 WRITE(numout,*) ' Mean surface error on big error points ', meance, numit 594 WRITE(numout,*) ' Number of points where there is a surf err gt than surf_err : ', numce, numit 595 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) 599 600 !--------------------------------------- 601 ! Write ice state in case of big errors 602 !--------------------------------------- 603 DO ji = kideb, kiut 604 IF ( ( ( t_su_b(ji) .LT. rtt ) .AND. ( surf_error(ji,jl) .GT. max_surf_err ) ) .OR. & 605 ( cons_error(ji,jl) .GT. max_cons_err ) ) THEN 606 ii = MOD( npb(ji) - 1, jpi ) + 1 607 ij = ( npb(ji) - 1 ) / jpi + 1 608 ! 609 WRITE(numout,*) ' alerte 1 ' 610 WRITE(numout,*) ' Untolerated conservation / surface error after ' 611 WRITE(numout,*) ' heat diffusion in the ice ' 612 WRITE(numout,*) ' Category : ', jl 613 WRITE(numout,*) ' ii , ij : ', ii, ij 614 WRITE(numout,*) ' lat, lon : ', gphit(ii,ij), glamt(ii,ij) 615 WRITE(numout,*) ' cons_error : ', cons_error(ji,jl) 616 WRITE(numout,*) ' surf_error : ', surf_error(ji,jl) 617 WRITE(numout,*) ' dq_i : ', - dq_i(ji,jl) * r1_rdtice 618 WRITE(numout,*) ' Fdt : ', sum_fluxq(ji,jl) 619 WRITE(numout,*) 620 ! WRITE(numout,*) ' qt_i_in : ', qt_i_in(ji,jl) 621 ! WRITE(numout,*) ' qt_s_in : ', qt_s_in(ji,jl) 622 ! WRITE(numout,*) ' qt_i_fin : ', qt_i_fin(ji,jl) 623 ! WRITE(numout,*) ' qt_s_fin : ', qt_s_fin(ji,jl) 624 ! WRITE(numout,*) ' qt : ', qt_i_fin(ji,jl) + qt_s_fin(ji,jl) 625 WRITE(numout,*) ' ht_i : ', ht_i_b(ji) 626 WRITE(numout,*) ' ht_s : ', ht_s_b(ji) 627 WRITE(numout,*) ' t_su : ', t_su_b(ji) 628 WRITE(numout,*) ' t_s : ', t_s_b(ji,1) 629 WRITE(numout,*) ' t_i : ', t_i_b(ji,1:nlay_i) 630 WRITE(numout,*) ' t_bo : ', t_bo_b(ji) 631 WRITE(numout,*) ' q_i : ', q_i_b(ji,1:nlay_i) 632 WRITE(numout,*) ' s_i : ', s_i_b(ji,1:nlay_i) 633 WRITE(numout,*) ' tmelts : ', rtt - tmut*s_i_b(ji,1:nlay_i) 634 WRITE(numout,*) 635 WRITE(numout,*) ' Fluxes ' 636 WRITE(numout,*) ' ~~~~~~ ' 637 WRITE(numout,*) ' fatm : ', fatm(ji,jl) 638 WRITE(numout,*) ' fc_su : ', fc_su (ji) 639 WRITE(numout,*) ' fstr_inice : ', qsr_ice_1d(ji)*i0(ji) 640 WRITE(numout,*) ' fc_bo : ', - fc_bo_i (ji) 641 WRITE(numout,*) ' foc : ', fbif_1d(ji) 642 WRITE(numout,*) ' fstroc : ', fstroc (ii,ij,jl) 643 WRITE(numout,*) ' i0 : ', i0(ji) 644 WRITE(numout,*) ' qsr_ice : ', (1.0-i0(ji))*qsr_ice_1d(ji) 645 WRITE(numout,*) ' qns_ice : ', qnsr_ice_1d(ji) 646 WRITE(numout,*) ' Conduction fluxes : ' 647 WRITE(numout,*) ' fc_s : ', fc_s(ji,0:nlay_s) 648 WRITE(numout,*) ' fc_i : ', fc_i(ji,0:nlay_i) 649 WRITE(numout,*) 650 WRITE(numout,*) ' Layer by layer ... ' 651 WRITE(numout,*) ' dq_snow : ', ( qt_s_fin(ji,jl) - qt_s_in(ji,jl) ) * r1_rdtice 652 WRITE(numout,*) ' dfc_snow : ', fc_s(ji,1) - fc_s(ji,0) 653 DO jk = 1, nlay_i 654 WRITE(numout,*) ' layer : ', jk 655 WRITE(numout,*) ' dq_ice : ', dq_i_layer(ji,jk) * r1_rdtice 656 WRITE(numout,*) ' radab : ', radab(ji,jk) 657 WRITE(numout,*) ' dfc_i : ', fc_i(ji,jk) - fc_i(ji,jk-1) 658 WRITE(numout,*) ' tot f : ', fc_i(ji,jk) - fc_i(ji,jk-1) - radab(ji,jk) 659 END DO 660 661 ENDIF 662 ! 663 END DO 664 ! 665 END SUBROUTINE lim_thd_con_dif 666 667 668 SUBROUTINE lim_thd_con_dh( kideb, kiut, jl ) 669 !!----------------------------------------------------------------------- 670 !! *** ROUTINE lim_thd_con_dh *** 671 !! 672 !! ** Purpose : Test energy conservation after enthalpy redistr. 673 !!----------------------------------------------------------------------- 674 INTEGER, INTENT(in) :: kideb, kiut ! bounds for the spatial loop 675 INTEGER, INTENT(in) :: jl ! category number 676 ! 677 INTEGER :: ji ! loop indices 678 INTEGER :: ii, ij, numce ! local integers 679 REAL(wp) :: meance, max_cons_err !local scalar 680 !!--------------------------------------------------------------------- 681 682 max_cons_err = 1._wp 683 684 !-------------------------- 685 ! Increment of energy 686 !-------------------------- 687 DO ji = kideb, kiut 688 dq_i(ji,jl) = qt_i_fin(ji,jl) - qt_i_in(ji,jl) + qt_s_fin(ji,jl) - qt_s_in(ji,jl) ! global 689 END DO 690 dq_i_layer(:,:) = q_i_layer_fin(:,:) - q_i_layer_in(:,:) ! layer by layer 691 692 !---------------------------------------- 693 ! Atmospheric heat flux, ice heat budget 694 !---------------------------------------- 695 DO ji = kideb, kiut 696 ii = MOD( npb(ji) - 1 , jpi ) + 1 697 ij = ( npb(ji) - 1 ) / jpi + 1 698 699 fatm (ji,jl) = qnsr_ice_1d(ji) + qsr_ice_1d(ji) ! total heat flux 700 sum_fluxq (ji,jl) = fatm(ji,jl) + fbif_1d(ji) - ftotal_fin(ji) - fstroc(ii,ij,jl) 701 cons_error(ji,jl) = ABS( dq_i(ji,jl) * r1_rdtice + sum_fluxq(ji,jl) ) 702 END DO 703 704 !-------------------- 705 ! Conservation error 706 !-------------------- 707 DO ji = kideb, kiut 708 cons_error(ji,jl) = ABS( dq_i(ji,jl) * r1_rdtice + sum_fluxq(ji,jl) ) 709 END DO 710 711 numce = 0 712 meance = 0._wp 713 DO ji = kideb, kiut 714 IF( cons_error(ji,jl) .GT. max_cons_err ) THEN 715 numce = numce + 1 716 meance = meance + cons_error(ji,jl) 717 ENDIF 718 ENDDO 719 IF(numce > 0 ) meance = meance / numce 720 721 WRITE(numout,*) ' Error report - Category : ', jl 722 WRITE(numout,*) ' ~~~~~~~~~~~~ ' 723 WRITE(numout,*) ' Maximum tolerated conservation error : ', max_cons_err 724 WRITE(numout,*) ' After lim_thd_ent, category : ', jl 725 WRITE(numout,*) ' Mean conservation error on big error points ', meance, numit 726 WRITE(numout,*) ' Number of points where there is a cons err gt than 0.1 W/m2 : ', numce, numit 727 728 !--------------------------------------- 729 ! Write ice state in case of big errors 730 !--------------------------------------- 731 DO ji = kideb, kiut 732 IF ( cons_error(ji,jl) .GT. max_cons_err ) THEN 733 ii = MOD( npb(ji) - 1, jpi ) + 1 734 ij = ( npb(ji) - 1 ) / jpi + 1 735 ! 736 WRITE(numout,*) ' alerte 1 - category : ', jl 737 WRITE(numout,*) ' Untolerated conservation error after limthd_ent ' 738 WRITE(numout,*) ' ii , ij : ', ii, ij 739 WRITE(numout,*) ' lat, lon : ', gphit(ii,ij), glamt(ii,ij) 740 WRITE(numout,*) ' * ' 741 WRITE(numout,*) ' Ftotal : ', sum_fluxq(ji,jl) 742 WRITE(numout,*) ' dq_t : ', - dq_i(ji,jl) * r1_rdtice 743 WRITE(numout,*) ' dq_i : ', - ( qt_i_fin(ji,jl) - qt_i_in(ji,jl) ) * r1_rdtice 744 WRITE(numout,*) ' dq_s : ', - ( qt_s_fin(ji,jl) - qt_s_in(ji,jl) ) * r1_rdtice 745 WRITE(numout,*) ' cons_error : ', cons_error(ji,jl) 746 WRITE(numout,*) ' * ' 747 WRITE(numout,*) ' Fluxes --- : ' 748 WRITE(numout,*) ' fatm : ', fatm(ji,jl) 749 WRITE(numout,*) ' foce : ', fbif_1d(ji) 750 WRITE(numout,*) ' fres : ', ftotal_fin(ji) 751 WRITE(numout,*) ' fhbri : ', fhbricat(ii,ij,jl) 752 WRITE(numout,*) ' * ' 753 WRITE(numout,*) ' Heat contents --- : ' 754 WRITE(numout,*) ' qt_s_in : ', qt_s_in(ji,jl) * r1_rdtice 755 WRITE(numout,*) ' qt_i_in : ', qt_i_in(ji,jl) * r1_rdtice 756 WRITE(numout,*) ' qt_in : ', ( qt_i_in(ji,jl) + qt_s_in(ji,jl) ) * r1_rdtice 757 WRITE(numout,*) ' qt_s_fin : ', qt_s_fin(ji,jl) * r1_rdtice 758 WRITE(numout,*) ' qt_i_fin : ', qt_i_fin(ji,jl) * r1_rdtice 759 WRITE(numout,*) ' qt_fin : ', ( qt_i_fin(ji,jl) + qt_s_fin(ji,jl) ) * r1_rdtice 760 WRITE(numout,*) ' * ' 761 WRITE(numout,*) ' Ice variables --- : ' 762 WRITE(numout,*) ' ht_i : ', ht_i_b(ji) 763 WRITE(numout,*) ' ht_s : ', ht_s_b(ji) 764 WRITE(numout,*) ' dh_s_tot : ', dh_s_tot(ji) 765 WRITE(numout,*) ' dh_snowice: ', dh_snowice(ji) 766 WRITE(numout,*) ' dh_i_surf : ', dh_i_surf(ji) 767 WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji) 768 ENDIF 769 ! 770 END DO 771 ! 772 END SUBROUTINE lim_thd_con_dh 773 774 555 775 556 SUBROUTINE lim_thd_enmelt( kideb, kiut ) 776 557 !!----------------------------------------------------------------------- … … 859 640 WRITE(numout,*)' maximal err. on T for heat diffusion computation maxer_i_thd = ', maxer_i_thd 860 641 WRITE(numout,*)' switch for comp. of thermal conductivity in the ice thcon_i_swi = ', thcon_i_swi 642 WRITE(numout,*)' check heat conservation in the ice/snow con_i = ', con_i 861 643 ENDIF 862 644 !
Note: See TracChangeset
for help on using the changeset viewer.
