Changeset 5965 for branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90
- Timestamp:
- 2015-12-01T16:35:30+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90
- Property svn:keywords set to Id
r4333 r5965 5 5 !!====================================================================== 6 6 !! History : 3.0 ! 2006-04 (M. Vancoppenolle) Original code 7 !! 3.5 ! 2014-06 (C. Rousset) Complete rewriting/cleaning 7 8 !!---------------------------------------------------------------------- 8 9 #if defined key_lim3 … … 12 13 !! lim_update1 : computes update of sea-ice global variables from trend terms 13 14 !!---------------------------------------------------------------------- 14 USE limrhg ! ice rheology15 16 USE dom_oce17 USE oce ! dynamics and tracers variables18 USE in_out_manager19 15 USE sbc_oce ! Surface boundary condition: ocean fields 20 16 USE sbc_ice ! Surface boundary condition: ice fields 21 17 USE dom_ice 18 USE dom_oce 22 19 USE phycst ! physical constants 23 20 USE ice 24 USE limdyn25 USE limtrp26 USE limthd27 USE limsbc28 USE limdiahsb29 USE limwri30 USE limrst31 21 USE thd_ice ! LIM thermodynamic sea-ice variables 32 USE par_ice33 22 USE limitd_th 34 23 USE limvar 35 USE prtctl ! Print control 36 USE lbclnk ! lateral boundary condition - MPP exchanges 37 USE wrk_nemo ! work arrays 38 USE lib_fortran ! glob_sum 39 ! Check budget (Rousset) 40 USE in_out_manager ! I/O manager 41 USE iom ! I/O manager 42 USE lib_mpp ! MPP library 24 USE prtctl ! Print control 25 USE wrk_nemo ! work arrays 43 26 USE timing ! Timing 27 USE limcons ! conservation tests 28 USE lib_mpp ! MPP library 29 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 30 USE in_out_manager ! I/O manager 44 31 45 32 IMPLICIT NONE 46 33 PRIVATE 47 34 48 PUBLIC lim_update1 ! routine called by ice_step 49 50 REAL(wp) :: epsi10 = 1.e-10_wp ! - - 51 REAL(wp) :: rzero = 0._wp ! - - 52 REAL(wp) :: rone = 1._wp ! - - 53 35 PUBLIC lim_update1 36 54 37 !! * Substitutions 55 38 # include "vectopt_loop_substitute.h90" 56 39 !!---------------------------------------------------------------------- 57 40 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 58 !! $Id : limupdate.F90 3294 2012-01-28 16:44:18Z rblod$41 !! $Id$ 59 42 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 60 43 !!---------------------------------------------------------------------- 61 44 CONTAINS 62 45 63 SUBROUTINE lim_update1 46 SUBROUTINE lim_update1( kt ) 64 47 !!------------------------------------------------------------------- 65 48 !! *** ROUTINE lim_update1 *** 66 49 !! 67 50 !! ** Purpose : Computes update of sea-ice global variables at 68 !! the end of the time step. 69 !! Address pathological cases 70 !! This place is very important 51 !! the end of the dynamics. 71 52 !! 72 !! ** Method :73 !! Ice speed from ice dynamics74 !! Ice thickness, Snow thickness, Temperatures, Lead fraction75 !! from advection and ice thermodynamics76 !!77 !! ** Action : -78 53 !!--------------------------------------------------------------------- 79 INTEGER :: ji, jj, jk, jl, jm ! dummy loop indices 80 INTEGER :: jbnd1, jbnd2 81 INTEGER :: i_ice_switch 82 INTEGER :: ind_im, layer ! indices for internal melt 83 REAL(wp) :: zweight, zesum, z_da_i, zhimax 84 REAL(wp) :: zinda, zindb, zindsn, zindic 85 REAL(wp) :: zindg, zh, zdvres, zviold2 86 REAL(wp) :: zbigvalue, zvsold2, z_da_ex 87 REAL(wp) :: z_prescr_hi, zat_i_old, ztmelts, ze_s 88 89 REAL(wp), POINTER, DIMENSION(:) :: zthick0, zqm0 ! thickness of the layers and heat contents for 90 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) 91 REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 92 ! mass and salt flux (clem) 93 REAL(wp), POINTER, DIMENSION(:,:,:) :: zviold, zvsold, zsmvold ! old ice volume... 54 INTEGER, INTENT(in) :: kt ! number of iteration 55 INTEGER :: ji, jj, jk, jl ! dummy loop indices 56 REAL(wp) :: zsal 57 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 94 58 !!------------------------------------------------------------------- 95 59 IF( nn_timing == 1 ) CALL timing_start('limupdate1') 96 60 97 CALL wrk_alloc( jkmax, zthick0, zqm0 ) 98 99 CALL wrk_alloc( jpi,jpj,jpl,zviold, zvsold, zsmvold ) ! clem 100 101 !------------------------------------------------------------------------------ 102 ! 1. Update of Global variables | 103 !------------------------------------------------------------------------------ 104 105 !----------------- 106 ! Trend terms 107 !----------------- 108 d_u_ice_dyn(:,:) = u_ice(:,:) - old_u_ice(:,:) 109 d_v_ice_dyn(:,:) = v_ice(:,:) - old_v_ice(:,:) 110 d_a_i_trp (:,:,:) = a_i (:,:,:) - old_a_i (:,:,:) 111 d_v_s_trp (:,:,:) = v_s (:,:,:) - old_v_s (:,:,:) 112 d_v_i_trp (:,:,:) = v_i (:,:,:) - old_v_i (:,:,:) 113 d_e_s_trp (:,:,:,:) = e_s (:,:,:,:) - old_e_s (:,:,:,:) 114 d_e_i_trp (:,:,:,:) = e_i (:,:,:,:) - old_e_i (:,:,:,:) 115 d_oa_i_trp (:,:,:) = oa_i (:,:,:) - old_oa_i (:,:,:) 116 d_smv_i_trp(:,:,:) = 0._wp 117 IF( num_sal == 2 ) d_smv_i_trp(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:) 118 119 ! mass and salt flux init (clem) 120 zviold(:,:,:) = v_i(:,:,:) 121 zvsold(:,:,:) = v_s(:,:,:) 122 zsmvold(:,:,:) = smv_i(:,:,:) 123 124 ! ------------------------------- 125 !- check conservation (C Rousset) 126 IF (ln_limdiahsb) THEN 127 zchk_v_i_b = glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 128 zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 129 zchk_fw_b = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) 130 zchk_fs_b = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) 61 IF( ln_limdyn ) THEN 62 63 IF( kt == nit000 .AND. lwp ) THEN 64 WRITE(numout,*) ' lim_update1 ' 65 WRITE(numout,*) ' ~~~~~~~~~~~ ' 131 66 ENDIF 132 !- check conservation (C Rousset) 133 ! ------------------------------- 134 135 CALL lim_var_glo2eqv 136 137 !-------------------------------------- 138 ! 2. Review of all pathological cases 139 !-------------------------------------- 140 141 ! clem: useless now 142 !------------------------------------------- 143 ! 2.1) Advection of ice in an ice-free cell 144 !------------------------------------------- 145 ! should be removed since it is treated after dynamics now 146 ! zhimax = 5._wp 147 ! ! first category 148 ! DO jj = 1, jpj 149 ! DO ji = 1, jpi 150 ! !--- the thickness of such an ice is often out of bounds 151 ! !--- thus we recompute a new area while conserving ice volume 152 ! zat_i_old = SUM( old_a_i(ji,jj,:) ) 153 ! zindb = MAX( 0._wp, SIGN( 1._wp, ABS( d_a_i_trp(ji,jj,1) ) - epsi10 ) ) 154 ! IF( ( ABS( d_v_i_trp(ji,jj,1) ) / MAX( ABS( d_a_i_trp(ji,jj,1) ), epsi10 ) * zindb .GT. zhimax ) & 155 ! & .AND.( ( v_i(ji,jj,1) / MAX( a_i(ji,jj,1), epsi10 ) * zindb ) .GT. zhimax ) & 156 ! & .AND.( zat_i_old .LT. 1.e-6 ) ) THEN ! new line 157 ! ht_i(ji,jj,1) = hi_max(1) * 0.5_wp 158 ! a_i (ji,jj,1) = v_i(ji,jj,1) / ht_i(ji,jj,1) 159 ! ENDIF 160 ! END DO 161 ! END DO 162 ! 163 ! zhimax = 20._wp 164 ! ! other categories 165 ! DO jl = 2, jpl 166 ! jm = ice_types(jl) 167 ! DO jj = 1, jpj 168 ! DO ji = 1, jpi 169 ! zindb = MAX( rzero, SIGN( rone, ABS( d_a_i_trp(ji,jj,jl) ) - epsi10 ) ) 170 ! ! this correction is very tricky... sometimes, advection gets wrong i don't know why 171 ! ! it makes problems when the advected volume and concentration do not seem to be 172 ! ! related with each other 173 ! ! the new thickness is sometimes very big! 174 ! ! and sometimes d_a_i_trp and d_v_i_trp have different sign 175 ! ! which of course is plausible 176 ! ! but fuck! it fucks everything up :) 177 ! IF ( ( ABS( d_v_i_trp(ji,jj,jl) ) / MAX( ABS( d_a_i_trp(ji,jj,jl) ), epsi10 ) * zindb .GT. zhimax ) & 178 ! & .AND. ( v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zindb ) .GT. zhimax ) THEN 179 ! ht_i(ji,jj,jl) = ( hi_max_typ(jl-ice_cat_bounds(jm,1),jm) + hi_max_typ(jl-ice_cat_bounds(jm,1)+1,jm) ) * 0.5_wp 180 ! a_i (ji,jj,jl) = v_i(ji,jj,jl) / ht_i(ji,jj,jl) 181 ! ENDIF 182 ! END DO ! ji 183 ! END DO !jj 184 ! END DO !jl 185 67 68 ! conservation test 69 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate1', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 70 71 !---------------------------------------------------- 72 ! ice concentration should not exceed amax 73 !----------------------------------------------------- 186 74 at_i(:,:) = 0._wp 187 75 DO jl = 1, jpl … … 189 77 END DO 190 78 191 !---------------------------------------------------- 192 ! 2.2) Rebin categories with thickness out of bounds 193 !---------------------------------------------------- 194 DO jm = 1, jpm 195 jbnd1 = ice_cat_bounds(jm,1) 196 jbnd2 = ice_cat_bounds(jm,2) 197 IF (ice_ncat_types(jm) .GT. 1 ) CALL lim_itd_th_reb(jbnd1, jbnd2, jm) 198 END DO 199 200 at_i(:,:) = 0._wp 201 DO jl = 1, jpl 202 at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 203 END DO 204 205 zbigvalue = 1.0e+20 206 207 DO jl = 1, jpl 208 DO jj = 1, jpj 79 DO jl = 1, jpl 80 DO jj = 1, jpj 209 81 DO ji = 1, jpi 210 211 !switches 212 zindb = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi10 ) ) 213 !switch = 1 if a_i > 1e-06 and 0 if not 214 zindsn = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) ) !=1 if hs > 1e-10 and 0 if not 215 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) !=1 if hi > 1e-10 and 0 if not 216 ! bug fix 25 avril 2007 217 zindb = zindb*zindic 218 219 !--- 2.3 Correction to ice age 220 !------------------------------ 221 ! IF ((o_i(ji,jj,jl)-1.0)*rday.gt.(rdt_ice*float(numit))) THEN 222 ! o_i(ji,jj,jl) = rdt_ice*FLOAT(numit)/rday 223 ! ENDIF 224 IF ((oa_i(ji,jj,jl)-1.0)*rday.gt.(rdt_ice*numit*a_i(ji,jj,jl))) THEN 225 oa_i(ji,jj,jl) = rdt_ice*numit/rday*a_i(ji,jj,jl) 82 IF( at_i(ji,jj) > rn_amax .AND. a_i(ji,jj,jl) > 0._wp ) THEN 83 a_i (ji,jj,jl) = a_i (ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax / at_i(ji,jj) ) ) 84 oa_i(ji,jj,jl) = oa_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax / at_i(ji,jj) ) ) 226 85 ENDIF 227 oa_i(ji,jj,jl) = zindb*zindic*oa_i(ji,jj,jl) 228 229 !--- 2.4 Correction to snow thickness 230 !------------------------------------- 231 ! ! snow thickness has to be greater than 0, and if ice concentration smaller than 1e-6 then hs = 0 232 ! v_s(ji,jj,jl) = MAX( zindb * v_s(ji,jj,jl), 0.0) 233 ! snow thickness cannot be smaller than 1e-6 234 zdvres = (zindsn * zindb - 1._wp) * v_s(ji,jj,jl) 235 v_s(ji,jj,jl) = v_s(ji,jj,jl) + zdvres 236 237 !rdm_snw(ji,jj) = rdm_snw(ji,jj) + zdvres * rhosn 238 239 !--- 2.5 Correction to ice thickness 240 !------------------------------------- 241 zdvres = (zindb - 1._wp) * v_i(ji,jj,jl) 242 v_i(ji,jj,jl) = v_i(ji,jj,jl) + zdvres 243 244 !rdm_ice(ji,jj) = rdm_ice(ji,jj) + zdvres * rhoic 245 !sfx_res(ji,jj) = sfx_res(ji,jj) - sm_i(ji,jj,jl) * ( rhoic * zdvres / rdt_ice ) 246 247 !--- 2.6 Snow is transformed into ice if the original ice cover disappears 248 !---------------------------------------------------------------------------- 249 zindg = tms(ji,jj) * MAX( 0._wp, SIGN( 1._wp, -v_i(ji,jj,jl) ) ) 250 zdvres = zindg * rhosn * v_s(ji,jj,jl) / rau0 251 v_i(ji,jj,jl) = v_i(ji,jj,jl) + zdvres 252 253 zdvres = zindsn*zindb * ( - zindg * v_s(ji,jj,jl) + zindg * v_i(ji,jj,jl) * ( rau0 - rhoic ) / rhosn ) 254 v_s(ji,jj,jl) = v_s(ji,jj,jl) + zdvres 255 256 !--- 2.7 Correction to ice concentrations 257 !-------------------------------------------- 258 ! if greater than 0, ice concentration cannot be smaller than 1e-10 259 a_i(ji,jj,jl) = zindb * a_i(ji,jj,jl) 260 261 !------------------------- 262 ! 2.8) Snow heat content 263 !------------------------- 264 e_s(ji,jj,1,jl) = zindsn * ( MIN ( MAX ( 0._wp, e_s(ji,jj,1,jl) ), zbigvalue ) ) 265 266 END DO ! ji 267 END DO ! jj 268 END DO ! jl 269 270 !------------------------ 271 ! 2.9) Ice heat content 272 !------------------------ 273 274 DO jl = 1, jpl 275 DO jk = 1, nlay_i 86 END DO 87 END DO 88 END DO 89 90 !--------------------- 91 ! Ice salinity bounds 92 !--------------------- 93 IF ( nn_icesal == 2 ) THEN 94 DO jl = 1, jpl 276 95 DO jj = 1, jpj 277 96 DO ji = 1, jpi 278 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) 279 e_i(ji,jj,jk,jl)= zindic * ( MIN ( MAX ( 0.0, e_i(ji,jj,jk,jl) ), zbigvalue ) ) 280 END DO ! ji 281 END DO ! jj 282 END DO !jk 283 END DO !jl 284 285 at_i(:,:) = 0._wp 286 DO jl = 1, jpl 287 at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 288 END DO 289 290 !--- 2.13 ice concentration should not exceed amax 291 ! (it should not be the case) 292 !----------------------------------------------------- 97 zsal = smv_i(ji,jj,jl) 98 ! salinity stays in bounds 99 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) ) 100 smv_i(ji,jj,jl) = rswitch * MAX( MIN( rn_simax * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), rn_simin * v_i(ji,jj,jl) ) 101 ! associated salt flux 102 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice 103 END DO 104 END DO 105 END DO 106 ENDIF 107 108 !---------------------------------------------------- 109 ! Rebin categories with thickness out of bounds 110 !---------------------------------------------------- 111 IF ( jpl > 1 ) CALL lim_itd_th_reb(1, jpl) 112 113 !----------------- 114 ! zap small values 115 !----------------- 116 CALL lim_var_zapsmall 117 118 ! ------------------------------------------------- 119 ! Diagnostics 120 ! ------------------------------------------------- 121 DO jl = 1, jpl 122 afx_dyn(:,:) = afx_dyn(:,:) + ( a_i(:,:,jl) - a_i_b(:,:,jl) ) * r1_rdtice 123 END DO 124 293 125 DO jj = 1, jpj 294 DO ji = 1, jpi 295 z_da_ex = MAX( at_i(ji,jj) - amax , 0.0 ) 296 zindb = MAX( rzero, SIGN( rone, at_i(ji,jj) - epsi10 ) ) 297 DO jl = 1, jpl 298 z_da_i = a_i(ji,jj,jl) * z_da_ex / MAX( at_i(ji,jj), epsi10 ) * zindb 299 a_i(ji,jj,jl) = MAX( 0._wp, a_i(ji,jj,jl) - z_da_i ) 300 ! 301 zinda = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi10 ) ) 302 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zinda 303 !v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) ! makes ice shrinken but should not be used 304 END DO 305 END DO 306 END DO 307 at_i(:,:) = a_i(:,:,1) 308 DO jl = 2, jpl 309 at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 310 END DO 311 312 313 ! Final thickness distribution rebinning 314 ! -------------------------------------- 315 DO jm = 1, jpm 316 jbnd1 = ice_cat_bounds(jm,1) 317 jbnd2 = ice_cat_bounds(jm,2) 318 IF (ice_ncat_types(jm) .GT. 1 ) CALL lim_itd_th_reb(jbnd1, jbnd2, jm) 319 IF (ice_ncat_types(jm) .EQ. 1 ) THEN 320 ENDIF 321 END DO 322 323 324 !--------------------- 325 ! 2.11) Ice salinity 326 !--------------------- 327 ! clem correct bug on smv_i 328 smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:) 329 330 IF ( num_sal == 2 ) THEN ! general case 331 DO jl = 1, jpl 332 !DO jk = 1, nlay_i 333 DO jj = 1, jpj 334 DO ji = 1, jpi 335 ! salinity stays in bounds 336 !clem smv_i(ji,jj,jl) = MAX(MIN((rhoic-rhosn)/rhoic*sss_m(ji,jj),smv_i(ji,jj,jl)),0.1 * v_i(ji,jj,jl) ) 337 smv_i(ji,jj,jl) = MAX( MIN( s_i_max * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), s_i_min * v_i(ji,jj,jl) ) 338 i_ice_switch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -v_i(ji,jj,jl) ) ) 339 smv_i(ji,jj,jl) = i_ice_switch * smv_i(ji,jj,jl) !+ s_i_min * ( 1._wp - i_ice_switch ) * v_i(ji,jj,jl) 340 END DO ! ji 341 END DO ! jj 342 !END DO !jk 343 END DO !jl 344 ENDIF 345 346 at_i(:,:) = a_i(:,:,1) 347 DO jl = 2, jpl 348 at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 349 END DO 350 351 352 !-------------------------------- 353 ! Update mass/salt fluxes (clem) 354 !-------------------------------- 355 DO jl = 1, jpl 356 DO jj = 1, jpj 357 DO ji = 1, jpi 358 diag_res_pr(ji,jj) = diag_res_pr(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) / rdt_ice 359 rdm_ice(ji,jj) = rdm_ice(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * rhoic 360 rdm_snw(ji,jj) = rdm_snw(ji,jj) + ( v_s(ji,jj,jl) - zvsold(ji,jj,jl) ) * rhosn 361 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmvold(ji,jj,jl) ) * rhoic / rdt_ice 362 END DO 363 END DO 364 END DO 365 366 ! ------------------------------- 367 !- check conservation (C Rousset) 368 IF (ln_limdiahsb) THEN 369 370 zchk_fs = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 371 zchk_fw = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b 372 373 zchk_v_i = ( glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) * r1_rdtice 374 zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic ) 375 376 zchk_vmin = glob_min(v_i) 377 zchk_amax = glob_max(SUM(a_i,dim=3)) 378 zchk_amin = glob_min(a_i) 379 380 IF(lwp) THEN 381 IF ( ABS( zchk_v_i ) > 1.e-5 ) WRITE(numout,*) 'violation volume [m3/day] (limupdate1) = ',(zchk_v_i * rday) 382 IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limupdate1) = ',(zchk_smv * rday) 383 IF ( zchk_vmin < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limupdate1) = ',(zchk_vmin * 1.e-3) 384 IF ( zchk_amax > amax+epsi10 ) WRITE(numout,*) 'violation a_i>amax (limupdate1) = ',zchk_amax 385 IF ( zchk_amin < 0. ) WRITE(numout,*) 'violation a_i<0 (limupdate1) = ',zchk_amin 386 ENDIF 387 ENDIF 388 !- check conservation (C Rousset) 389 ! ------------------------------- 390 126 DO ji = 1, jpi 127 ! heat content variation (W.m-2) 128 diag_heat(ji,jj) = - ( SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) ) + & 129 & SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) ) & 130 & ) * r1_rdtice 131 ! salt, volume 132 diag_smvi(ji,jj) = SUM( smv_i(ji,jj,:) - smv_i_b(ji,jj,:) ) * rhoic * r1_rdtice 133 diag_vice(ji,jj) = SUM( v_i (ji,jj,:) - v_i_b (ji,jj,:) ) * rhoic * r1_rdtice 134 diag_vsnw(ji,jj) = SUM( v_s (ji,jj,:) - v_s_b (ji,jj,:) ) * rhosn * r1_rdtice 135 END DO 136 END DO 137 138 ! conservation test 139 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limupdate1', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 140 141 ! ------------------------------------------------- 142 ! control prints 143 ! ------------------------------------------------- 391 144 IF(ln_ctl) THEN ! Control print 392 145 CALL prt_ctl_info(' ') 393 146 CALL prt_ctl_info(' - Cell values : ') 394 147 CALL prt_ctl_info(' ~~~~~~~~~~~~~ ') 395 CALL prt_ctl(tab2d_1= area, clinfo1=' lim_update1 : cell area :')148 CALL prt_ctl(tab2d_1=e12t , clinfo1=' lim_update1 : cell area :') 396 149 CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_update1 : at_i :') 397 150 CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_update1 : vt_i :') … … 399 152 CALL prt_ctl(tab2d_1=strength , clinfo1=' lim_update1 : strength :') 400 153 CALL prt_ctl(tab2d_1=u_ice , clinfo1=' lim_update1 : u_ice :', tab2d_2=v_ice , clinfo2=' v_ice :') 401 CALL prt_ctl(tab2d_1=d_u_ice_dyn, clinfo1=' lim_update1 : d_u_ice_dyn :', tab2d_2=d_v_ice_dyn, clinfo2=' d_v_ice_dyn :') 402 CALL prt_ctl(tab2d_1=old_u_ice , clinfo1=' lim_update1 : old_u_ice :', tab2d_2=old_v_ice , clinfo2=' old_v_ice :') 154 CALL prt_ctl(tab2d_1=u_ice_b , clinfo1=' lim_update1 : u_ice_b :', tab2d_2=v_ice_b , clinfo2=' v_ice_b :') 403 155 404 156 DO jl = 1, jpl … … 413 165 CALL prt_ctl(tab2d_1=o_i (:,:,jl) , clinfo1= ' lim_update1 : o_i : ') 414 166 CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_update1 : a_i : ') 415 CALL prt_ctl(tab2d_1=old_a_i (:,:,jl) , clinfo1= ' lim_update1 : old_a_i : ') 416 CALL prt_ctl(tab2d_1=d_a_i_trp (:,:,jl) , clinfo1= ' lim_update1 : d_a_i_trp : ') 167 CALL prt_ctl(tab2d_1=a_i_b (:,:,jl) , clinfo1= ' lim_update1 : a_i_b : ') 417 168 CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_update1 : v_i : ') 418 CALL prt_ctl(tab2d_1=old_v_i (:,:,jl) , clinfo1= ' lim_update1 : old_v_i : ') 419 CALL prt_ctl(tab2d_1=d_v_i_trp (:,:,jl) , clinfo1= ' lim_update1 : d_v_i_trp : ') 169 CALL prt_ctl(tab2d_1=v_i_b (:,:,jl) , clinfo1= ' lim_update1 : v_i_b : ') 420 170 CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_update1 : v_s : ') 421 CALL prt_ctl(tab2d_1=old_v_s (:,:,jl) , clinfo1= ' lim_update1 : old_v_s : ') 422 CALL prt_ctl(tab2d_1=d_v_s_trp (:,:,jl) , clinfo1= ' lim_update1 : d_v_s_trp : ') 423 CALL prt_ctl(tab2d_1=e_i (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1 : e_i1 : ') 424 CALL prt_ctl(tab2d_1=old_e_i (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1 : old_e_i1 : ') 425 CALL prt_ctl(tab2d_1=d_e_i_trp (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1 : de_i1_trp : ') 426 CALL prt_ctl(tab2d_1=e_i (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1 : e_i2 : ') 427 CALL prt_ctl(tab2d_1=old_e_i (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1 : old_e_i2 : ') 428 CALL prt_ctl(tab2d_1=d_e_i_trp (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1 : de_i2_trp : ') 171 CALL prt_ctl(tab2d_1=v_s_b (:,:,jl) , clinfo1= ' lim_update1 : v_s_b : ') 172 CALL prt_ctl(tab2d_1=e_i (:,:,1,jl) , clinfo1= ' lim_update1 : e_i1 : ') 173 CALL prt_ctl(tab2d_1=e_i_b (:,:,1,jl) , clinfo1= ' lim_update1 : e_i1_b : ') 174 CALL prt_ctl(tab2d_1=e_i (:,:,2,jl) , clinfo1= ' lim_update1 : e_i2 : ') 175 CALL prt_ctl(tab2d_1=e_i_b (:,:,2,jl) , clinfo1= ' lim_update1 : e_i2_b : ') 429 176 CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_update1 : e_snow : ') 430 CALL prt_ctl(tab2d_1=old_e_s (:,:,1,jl) , clinfo1= ' lim_update1 : old_e_snow : ') 431 CALL prt_ctl(tab2d_1=d_e_s_trp (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1 : d_e_s_trp : ') 177 CALL prt_ctl(tab2d_1=e_s_b (:,:,1,jl) , clinfo1= ' lim_update1 : e_snow_b : ') 432 178 CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_update1 : smv_i : ') 433 CALL prt_ctl(tab2d_1=old_smv_i (:,:,jl) , clinfo1= ' lim_update1 : old_smv_i : ') 434 CALL prt_ctl(tab2d_1=d_smv_i_trp(:,:,jl) , clinfo1= ' lim_update1 : d_smv_i_trp : ') 179 CALL prt_ctl(tab2d_1=smv_i_b (:,:,jl) , clinfo1= ' lim_update1 : smv_i_b : ') 435 180 CALL prt_ctl(tab2d_1=oa_i (:,:,jl) , clinfo1= ' lim_update1 : oa_i : ') 436 CALL prt_ctl(tab2d_1=old_oa_i (:,:,jl) , clinfo1= ' lim_update1 : old_oa_i : ') 437 CALL prt_ctl(tab2d_1=d_oa_i_trp (:,:,jl) , clinfo1= ' lim_update1 : d_oa_i_trp : ') 181 CALL prt_ctl(tab2d_1=oa_i_b (:,:,jl) , clinfo1= ' lim_update1 : oa_i_b : ') 438 182 439 183 DO jk = 1, nlay_i … … 446 190 CALL prt_ctl_info(' - Heat / FW fluxes : ') 447 191 CALL prt_ctl_info(' ~~~~~~~~~~~~~~~~~~ ') 448 CALL prt_ctl(tab2d_1=fmmec , clinfo1= ' lim_update1 : fmmec : ', tab2d_2=fhmec , clinfo2= ' fhmec : ')449 192 CALL prt_ctl(tab2d_1=sst_m , clinfo1= ' lim_update1 : sst : ', tab2d_2=sss_m , clinfo2= ' sss : ') 450 CALL prt_ctl(tab2d_1=fhbri , clinfo1= ' lim_update1 : fhbri : ', tab2d_2=fheat_mec , clinfo2= ' fheat_mec : ')451 193 452 194 CALL prt_ctl_info(' ') … … 458 200 ENDIF 459 201 460 461 CALL wrk_dealloc( jkmax, zthick0, zqm0 ) 462 463 CALL wrk_dealloc( jpi,jpj,jpl,zviold, zvsold, zsmvold ) ! clem 202 ENDIF ! ln_limdyn 464 203 465 204 IF( nn_timing == 1 ) CALL timing_stop('limupdate1')
Note: See TracChangeset
for help on using the changeset viewer.
