- Timestamp:
- 2014-05-27T11:28:12+02:00 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90
r4634 r4649 39 39 USE lib_fortran ! glob_sum 40 40 USE timing ! Timing 41 USE limcons ! conservation tests 41 42 42 43 IMPLICIT NONE … … 84 85 REAL(wp) :: zdE ! specific enthalpy difference (J/kg) 85 86 REAL(wp) :: zfmdt ! exchange mass flux x time step (J/m2), >0 towards the ocean 86 87 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 ! Check conservation (C Rousset) 88 REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 87 ! 88 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 89 89 !!------------------------------------------------------------------- 90 90 IF( nn_timing == 1 ) CALL timing_start('limupdate2') 91 91 92 !---------------------------------------------------------------------------------------- 93 ! 1. Computation of trend terms 94 !---------------------------------------------------------------------------------------- 95 96 ! ------------------------------- 97 !- check conservation (C Rousset) 98 IF (ln_limdiahsb) THEN 99 zchk_v_i_b = glob_sum( SUM( v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) 100 zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 101 zchk_e_i_b = glob_sum( SUM( e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) 102 zchk_fw_b = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) ) 103 zchk_fs_b = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) 104 zchk_ft_b = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) 105 ENDIF 106 !- check conservation (C Rousset) 107 ! ------------------------------- 108 92 ! conservation test 93 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate2', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 94 95 !----------------- 109 96 ! zap small values 110 97 !----------------- … … 113 100 CALL lim_var_glo2eqv 114 101 115 !--------------------------------------116 ! 2. Review of all pathological cases117 !--------------------------------------118 at_i(:,:) = 0._wp119 DO jl = 1, jpl120 at_i(:,:) = a_i(:,:,jl) + at_i(:,:)121 END DO122 123 102 !---------------------------------------------------- 124 ! 2.2)Rebin categories with thickness out of bounds103 ! Rebin categories with thickness out of bounds 125 104 !---------------------------------------------------- 126 105 DO jm = 1, jpm … … 130 109 END DO 131 110 132 133 !clem debug: it is done in limthd_dh now 134 ! ! Melt of snow 135 ! !-------------- 136 ! DO jl = 1, jpl 137 ! DO jj = 1, jpj 138 ! DO ji = 1, jpi 139 ! IF( v_s(ji,jj,jl) >= epsi20 ) THEN 140 ! ! If snow energy of melting smaller then Lf 141 ! ! Then all snow melts and heat go to the ocean 142 ! !IF ( zEs <= lfus ) THEN 143 ! IF( t_s(ji,jj,1,jl) >= rtt ) THEN 144 ! zdvres = - v_s(ji,jj,jl) 145 ! zEs = - e_s(ji,jj,1,jl) * unit_fac / ( area(ji,jj) * MAX( v_s(ji,jj,jl), epsi20 ) ) ! snow energy of melting (J.m-3) 146 ! ! Contribution to heat flux to the ocean [W.m-2], < 0 147 ! hfx_res(ji,jj) = hfx_res(ji,jj) - zEs * zdvres * r1_rdtice 148 ! ! Contribution to mass flux 149 ! wfx_snw(ji,jj) = wfx_snw(ji,jj) + rhosn * zdvres * r1_rdtice 150 ! ! updates 151 ! v_s (ji,jj,jl) = 0._wp 152 ! ht_s(ji,jj,jl) = 0._wp 153 ! e_s (ji,jj,1,jl) = 0._wp 154 ! t_s (ji,jj,1,jl) = rtt 155 ! ENDIF 156 ! ENDIF 157 ! END DO 158 ! END DO 159 ! END DO 160 !clem debug 161 162 !--- 2.12 Constrain the thickness of the smallest category above 10 cm 111 !---------------------------------------------------------------------- 112 ! Constrain the thickness of the smallest category above hiclim 163 113 !---------------------------------------------------------------------- 164 114 DO jm = 1, jpm … … 176 126 END DO !jm 177 127 178 !--- 2.13 ice concentration should not exceed amax179 128 !----------------------------------------------------- 180 at_i(:,:) = 0.0 129 ! ice concentration should not exceed amax 130 !----------------------------------------------------- 131 at_i(:,:) = 0._wp 181 132 DO jl = 1, jpl 182 133 at_i(:,:) = a_i(:,:,jl) + at_i(:,:) … … 199 150 END DO 200 151 152 ! -------------------------------------- 201 153 ! Final thickness distribution rebinning 202 154 ! -------------------------------------- … … 209 161 END DO 210 162 163 !----------------- 211 164 ! zap small values 212 165 !----------------- … … 232 185 ENDIF 233 186 234 235 ! -------------------236 at_i(:,:) = a_i(:,:,1)237 DO jl = 2, jpl238 at_i(:,:) = a_i(:,:,jl) + at_i(:,:)239 END DO240 241 187 !------------------------------------------------------------------------------ 242 188 ! 2) Corrections to avoid wrong values | … … 246 192 DO jj = 2, jpjm1 247 193 DO ji = 2, jpim1 248 IF ( at_i(ji,jj) .EQ. 0.0) THEN ! what to do if there is no ice249 IF ( at_i(ji+1,jj) .EQ. 0.0 ) u_ice(ji,jj) = 0.0! right side250 IF ( at_i(ji-1,jj) .EQ. 0.0 ) u_ice(ji-1,jj) = 0.0! left side251 IF ( at_i(ji,jj+1) .EQ. 0.0 ) v_ice(ji,jj) = 0.0! upper side252 IF ( at_i(ji,jj-1) .EQ. 0.0 ) v_ice(ji,jj-1) = 0.0! bottom side194 IF ( at_i(ji,jj) == 0._wp ) THEN ! what to do if there is no ice 195 IF ( at_i(ji+1,jj) == 0._wp ) u_ice(ji,jj) = 0._wp ! right side 196 IF ( at_i(ji-1,jj) == 0._wp ) u_ice(ji-1,jj) = 0._wp ! left side 197 IF ( at_i(ji,jj+1) == 0._wp ) v_ice(ji,jj) = 0._wp ! upper side 198 IF ( at_i(ji,jj-1) == 0._wp ) v_ice(ji,jj-1) = 0._wp ! bottom side 253 199 ENDIF 254 200 END DO … … 261 207 v_ice(:,:) = v_ice(:,:) * tmv(:,:) 262 208 263 264 209 ! ------------------------------------------------- 265 210 ! Diagnostics … … 276 221 dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday 277 222 278 ! ------------------------------- 279 !- check conservation (C Rousset) 280 IF( ln_limdiahsb ) THEN 281 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 282 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 283 zchk_ft = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) - zchk_ft_b 284 285 zchk_v_i = ( glob_sum( SUM( v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b ) * r1_rdtice - zchk_fw 286 zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic ) 287 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 288 289 zchk_vmin = glob_min(v_i) 290 zchk_amax = glob_max(SUM(a_i,dim=3)) 291 zchk_amin = glob_min(a_i) 292 293 IF(lwp) THEN 294 IF ( ABS( zchk_v_i ) > 1.e-4 ) WRITE(numout,*) 'violation volume [kg/day] (limupdate2) = ',(zchk_v_i * rday) 295 IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limupdate2) = ',(zchk_smv * rday) 296 IF ( ABS( zchk_e_i ) > 1.e-2 ) WRITE(numout,*) 'violation enthalpy [1e9 J] (limupdate2) = ',(zchk_e_i) 297 IF ( zchk_vmin < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limupdate2) = ',(zchk_vmin * 1.e-3) 298 IF ( zchk_amax > amax+epsi10 ) WRITE(numout,*) 'violation a_i>amax (limupdate2) = ',zchk_amax 299 IF ( zchk_amin < 0. ) WRITE(numout,*) 'violation a_i<0 (limupdate2) = ',zchk_amin 300 ENDIF 301 ENDIF 302 !- check conservation (C Rousset) 303 ! ------------------------------- 223 ! heat content variation (W.m-2) 224 DO jj = 1, jpj 225 DO ji = 1, jpi 226 diag_heat_dhc(ji,jj) = ( SUM( d_e_i_trp(ji,jj,1:nlay_i,:) + d_e_i_thd(ji,jj,1:nlay_i,:) ) + & 227 & SUM( d_e_s_trp(ji,jj,1:nlay_s,:) + d_e_s_thd(ji,jj,1:nlay_s,:) ) ) * unit_fac * r1_rdtice / area(ji,jj) 228 END DO 229 END DO 230 231 ! conservation test 232 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate2', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 304 233 305 234 IF(ln_ctl) THEN ! Control print … … 370 299 371 300 IF( nn_timing == 1 ) CALL timing_stop('limupdate2') 301 372 302 END SUBROUTINE lim_update2 373 303 #else
Note: See TracChangeset
for help on using the changeset viewer.
