Changeset 4688 for trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
- Timestamp:
- 2014-06-25T01:39:59+02:00 (10 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90 (modified) (19 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
r4333 r4688 35 35 USE lib_fortran ! to use key_nosignedzero 36 36 USE timing ! Timing 37 USE limcons ! conservation tests 37 38 38 39 IMPLICIT NONE … … 65 66 INTEGER, INTENT(in) :: kt ! time step index 66 67 ! 67 INTEGER :: j l, ja, jm, jbnd1, jbnd2 ! ice types dummy loop index68 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)69 REAL(wp) :: z chk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset)68 INTEGER :: ji,jj, jk, jl, ja, jm, jbnd1, jbnd2 ! ice types dummy loop index 69 ! 70 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 70 71 !!------------------------------------------------------------------ 71 72 IF( nn_timing == 1 ) CALL timing_start('limitd_th') 72 73 73 ! ------------------------------- 74 !- check conservation (C Rousset) 75 IF (ln_limdiahsb) THEN 76 zchk_v_i_b = glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 77 zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 78 zchk_fw_b = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) 79 zchk_fs_b = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) 80 ENDIF 81 !- check conservation (C Rousset) 82 ! ------------------------------- 74 ! conservation test 75 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_th', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 83 76 84 77 IF( kt == nit000 .AND. lwp ) THEN … … 105 98 ! 3) Add frazil ice growing in leads. 106 99 !------------------------------------------------------------------------------| 107 108 100 CALL lim_thd_lac 109 101 CALL lim_var_glo2eqv ! only for info 110 111 IF(ln_ctl) THEN ! Control print102 103 IF(ln_ctl) THEN ! Control print 112 104 CALL prt_ctl_info(' ') 113 105 CALL prt_ctl_info(' - Cell values : ') … … 141 133 ENDIF 142 134 ! 143 ! ------------------------------- 144 !- check conservation (C Rousset) 145 IF( ln_limdiahsb ) THEN 146 zchk_fs = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 147 zchk_fw = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b 148 149 zchk_v_i = ( glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) * r1_rdtice 150 zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic ) 151 152 zchk_vmin = glob_min(v_i) 153 zchk_amax = glob_max(SUM(a_i,dim=3)) 154 zchk_amin = glob_min(a_i) 155 156 IF(lwp) THEN 157 IF ( ABS( zchk_v_i ) > 1.e-5 ) WRITE(numout,*) 'violation volume [m3/day] (limitd_th) = ',(zchk_v_i * rday) 158 IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limitd_th) = ',(zchk_smv * rday) 159 IF ( zchk_vmin < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limitd_th) = ',(zchk_vmin * 1.e-3) 160 IF ( zchk_amax > amax+epsi10 ) WRITE(numout,*) 'violation a_i>amax (limitd_th) = ',zchk_amax 161 IF ( zchk_amin < 0. ) WRITE(numout,*) 'violation a_i<0 (limitd_th) = ',zchk_amin 162 ENDIF 163 ENDIF 164 !- check conservation (C Rousset) 165 ! ------------------------------- 135 ! conservation test 136 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_th', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 166 137 ! 167 138 IF( nn_timing == 1 ) CALL timing_stop('limitd_th') … … 258 229 zindb = 1.0 - MAX( 0.0, SIGN( 1.0, - old_a_i(ji,jj,jl) + epsi10) ) !0 if no ice and 1 if yes 259 230 zht_i_o(ji,jj,jl) = old_v_i(ji,jj,jl) / MAX( old_a_i(ji,jj,jl), epsi10 ) * zindb 260 IF( a_i(ji,jj,jl) > epsi 06) zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_o(ji,jj,jl)231 IF( a_i(ji,jj,jl) > epsi10 ) zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_o(ji,jj,jl) 261 232 END DO 262 233 END DO … … 302 273 ij = nind_j(ji) 303 274 ! 304 IF ( ( zht_i_o(ii,ij,jl) .GT. epsi10 ) .AND. &305 ( zht_i_o(ii,ij,jl+1) .GT. epsi10 )) THEN275 zhbnew(ii,ij,jl) = hi_max(jl) 276 IF ( old_a_i(ii,ij,jl) > epsi10 .AND. old_a_i(ii,ij,jl+1) > epsi10 ) THEN 306 277 !interpolate between adjacent category growth rates 307 zslope = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / & 308 ( zht_i_o (ii,ij,jl+1) - zht_i_o (ii,ij,jl) ) 309 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + & 310 zslope * ( hi_max(jl) - zht_i_o(ii,ij,jl) ) 311 ELSEIF (zht_i_o(ii,ij,jl).gt.epsi10) THEN 278 zslope = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / ( zht_i_o(ii,ij,jl+1) - zht_i_o(ii,ij,jl) ) 279 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + zslope * ( hi_max(jl) - zht_i_o(ii,ij,jl) ) 280 ELSEIF ( old_a_i(ii,ij,jl) > epsi10) THEN 312 281 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) 313 ELSEIF ( zht_i_o(ii,ij,jl+1).gt.epsi10) THEN282 ELSEIF ( old_a_i(ii,ij,jl+1) > epsi10) THEN 314 283 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1) 315 ELSE316 zhbnew(ii,ij,jl) = hi_max(jl)317 284 ENDIF 318 285 END DO … … 320 287 !- 4.2 Check that each zhbnew lies between adjacent values of ice thickness 321 288 DO ji = 1, nbrem 322 ! jl, ji323 289 ii = nind_i(ji) 324 290 ij = nind_j(ji) 325 ! jl, ji 326 IF ( ( a_i(ii,ij,jl) .GT.epsi10) .AND. & 327 ( ht_i(ii,ij,jl).GE. zhbnew(ii,ij,jl) ) & 328 ) THEN 291 IF( a_i(ii,ij,jl) > epsi10 .AND. ht_i(ii,ij,jl) >= zhbnew(ii,ij,jl) ) THEN 329 292 zremap_flag(ii,ij) = 0 330 ELSEIF ( ( a_i(ii,ij,jl+1) .GT. epsi10 ) .AND. & 331 ( ht_i(ii,ij,jl+1).LE. zhbnew(ii,ij,jl) ) & 332 ) THEN 293 ELSEIF( a_i(ii,ij,jl+1) > epsi10 .AND. ht_i(ii,ij,jl+1) <= zhbnew(ii,ij,jl) ) THEN 333 294 zremap_flag(ii,ij) = 0 334 295 ENDIF 335 296 336 297 !- 4.3 Check that each zhbnew does not exceed maximal values hi_max 337 ! jl, ji 338 IF (zhbnew(ii,ij,jl).gt.hi_max(jl+1)) THEN 339 zremap_flag(ii,ij) = 0 340 ENDIF 341 ! jl, ji 342 IF (zhbnew(ii,ij,jl).lt.hi_max(jl-1)) THEN 343 zremap_flag(ii,ij) = 0 344 ENDIF 345 ! jl, ji 346 END DO !ji 347 ! ji 298 IF( zhbnew(ii,ij,jl) > hi_max(jl+1) ) zremap_flag(ii,ij) = 0 299 IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0 300 END DO 301 348 302 END DO !jl 349 303 … … 354 308 DO jj = 1, jpj 355 309 DO ji = 1, jpi 356 IF 310 IF( zremap_flag(ji,jj) == 1 ) THEN 357 311 nbrem = nbrem + 1 358 312 nind_i(nbrem) = ji 359 313 nind_j(nbrem) = jj 360 314 ENDIF 361 END DO !ji362 END DO !jj315 END DO 316 END DO 363 317 364 318 !----------------------------------------------------------------------------------------------- … … 380 334 ENDIF 381 335 382 IF( zhbnew(ji,jj,kubnd) < hi_max(kubnd-1) ) 336 IF( zhbnew(ji,jj,kubnd) < hi_max(kubnd-1) ) zhbnew(ji,jj,kubnd) = hi_max(kubnd-1) 383 337 384 338 END DO !jj … … 444 398 DO jl = klbnd, kubnd 445 399 CALL lim_itd_fitline(jl, zhbnew(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl), & 446 g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), & 447 zremap_flag) 400 g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), zremap_flag) 448 401 END DO 449 402 … … 493 446 nd = zdonor(ii,ij,jl) 494 447 zdaice(ii,ij,jl) = g1(ii,ij,nd)*zx2 + g0(ii,ij,nd)*zx1 495 zdvice(ii,ij,jl) = g1(ii,ij,nd)*zx3 + g0(ii,ij,nd)*zx2 + & 496 zdaice(ii,ij,jl)*hL(ii,ij,nd) 448 zdvice(ii,ij,jl) = g1(ii,ij,nd)*zx3 + g0(ii,ij,nd)*zx2 + zdaice(ii,ij,jl)*hL(ii,ij,nd) 497 449 498 450 END DO ! ji … … 511 463 ii = nind_i(ji) 512 464 ij = nind_j(ji) 513 IF ( ( a_i(ii,ij,1) > epsi10 ) .AND. ( ht_i(ii,ij,1) < hiclim )) THEN465 IF ( a_i(ii,ij,1) > epsi10 .AND. ht_i(ii,ij,1) < hiclim ) THEN 514 466 a_i(ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) / hiclim 515 467 ht_i(ii,ij,1) = hiclim 516 v_i(ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) !clem-useless517 468 ENDIF 518 469 END DO !ji … … 799 750 !-------------- 800 751 801 zdvsnow 752 zdvsnow = v_s(ii,ij,jl1) * zworka(ii,ij) 802 753 v_s(ii,ij,jl1) = v_s(ii,ij,jl1) - zdvsnow 803 754 v_s(ii,ij,jl2) = v_s(ii,ij,jl2) + zdvsnow … … 807 758 !-------------------- 808 759 809 zdesnow 760 zdesnow = e_s(ii,ij,1,jl1) * zworka(ii,ij) 810 761 e_s(ii,ij,1,jl1) = e_s(ii,ij,1,jl1) - zdesnow 811 762 e_s(ii,ij,1,jl2) = e_s(ii,ij,1,jl2) + zdesnow … … 815 766 !-------------- 816 767 817 zdo_aice 768 zdo_aice = oa_i(ii,ij,jl1) * zdaice(ii,ij,jl) 818 769 oa_i(ii,ij,jl1) = oa_i(ii,ij,jl1) - zdo_aice 819 770 oa_i(ii,ij,jl2) = oa_i(ii,ij,jl2) + zdo_aice … … 823 774 !-------------- 824 775 825 zdsm_vice 776 zdsm_vice = smv_i(ii,ij,jl1) * zworka(ii,ij) 826 777 smv_i(ii,ij,jl1) = smv_i(ii,ij,jl1) - zdsm_vice 827 778 smv_i(ii,ij,jl2) = smv_i(ii,ij,jl2) + zdsm_vice … … 831 782 !--------------------- 832 783 833 zdaTsf 784 zdaTsf = t_su(ii,ij,jl1) * zdaice(ii,ij,jl) 834 785 zaTsfn(ii,ij,jl1) = zaTsfn(ii,ij,jl1) - zdaTsf 835 786 zaTsfn(ii,ij,jl2) = zaTsfn(ii,ij,jl2) + zdaTsf … … 910 861 REAL(wp), POINTER, DIMENSION(:,:) :: vt_s_init, vt_s_final ! snow volume summed over categories 911 862 !!------------------------------------------------------------------ 863 !! clem 2014/04: be carefull, rebining does not conserve salt(maybe?) => the difference is taken into account in limupdate 912 864 913 865 CALL wrk_alloc( jpi,jpj,jpl, zdonor ) ! interger … … 1015 967 1016 968 !clem-change 969 DO jj = 1, jpj 970 DO ji = 1, jpi 971 IF( a_i(ji,jj,jl+1) > epsi10 .AND. ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN 972 ! 973 zshiftflag = 1 974 zdonor(ji,jj,jl) = jl + 1 975 zdaice(ji,jj,jl) = a_i(ji,jj,jl+1) 976 zdvice(ji,jj,jl) = v_i(ji,jj,jl+1) 977 ENDIF 978 END DO ! ji 979 END DO ! jj 980 981 IF(lk_mpp) CALL mpp_max( zshiftflag ) 982 983 IF( zshiftflag == 1 ) THEN ! Shift ice between categories 984 CALL lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice ) 985 ! Reset shift parameters 986 zdonor(:,:,jl) = 0 987 zdaice(:,:,jl) = 0._wp 988 zdvice(:,:,jl) = 0._wp 989 ENDIF 990 !clem-change 991 992 ! ! clem-change begin: why not doing that? 1017 993 ! DO jj = 1, jpj 1018 994 ! DO ji = 1, jpi 1019 ! IF( a_i(ji,jj,jl+1) > epsi10 .AND. & 1020 ! ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN 1021 ! ! 1022 ! zshiftflag = 1 1023 ! zdonor(ji,jj,jl) = jl + 1 1024 ! zdaice(ji,jj,jl) = a_i(ji,jj,jl+1) 1025 ! zdvice(ji,jj,jl) = v_i(ji,jj,jl+1) 995 ! IF( a_i(ji,jj,jl+1) > epsi10 .AND. ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN 996 ! ht_i(ji,jj,jl+1) = hi_max(jl) + epsi10 997 ! a_i (ji,jj,jl+1) = v_i(ji,jj,jl+1) / ht_i(ji,jj,jl+1) 1026 998 ! ENDIF 1027 999 ! END DO ! ji 1028 1000 ! END DO ! jj 1029 !1030 ! IF(lk_mpp) CALL mpp_max( zshiftflag )1031 !1032 ! IF( zshiftflag == 1 ) THEN ! Shift ice between categories1033 ! CALL lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice )1034 ! ! Reset shift parameters1035 ! zdonor(:,:,jl) = 01036 ! zdaice(:,:,jl) = 0._wp1037 ! zdvice(:,:,jl) = 0._wp1038 ! ENDIF1039 !clem-change1040 1041 ! clem-change begin: why not doing that?1042 DO jj = 1, jpj1043 DO ji = 1, jpi1044 IF( a_i(ji,jj,jl+1) > epsi10 .AND. &1045 ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN1046 ht_i(ji,jj,jl+1) = hi_max(jl) + epsi101047 a_i (ji,jj,jl+1) = v_i(ji,jj,jl+1) / ht_i(ji,jj,jl+1)1048 ENDIF1049 END DO ! ji1050 END DO ! jj1051 1001 ! clem-change end 1052 1002
Note: See TracChangeset
for help on using the changeset viewer.
