- 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/limitd_th.F90
r4332 r4634 65 65 INTEGER, INTENT(in) :: kt ! time step index 66 66 ! 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)67 INTEGER :: ji,jj, jk, jl, ja, jm, jbnd1, jbnd2 ! ice types dummy loop index 68 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) 69 69 REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 70 70 !!------------------------------------------------------------------ … … 74 74 !- check conservation (C Rousset) 75 75 IF (ln_limdiahsb) THEN 76 zchk_v_i_b = glob_sum( SUM( v_i(:,:,:) , dim=3 ) * area(:,:) * tms(:,:) )76 zchk_v_i_b = glob_sum( SUM( v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) 77 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(:,:) ) 78 zchk_e_i_b = glob_sum( SUM( e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) 79 zchk_fw_b = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) ) 80 zchk_fs_b = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) 81 zchk_ft_b = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) 80 82 ENDIF 81 83 !- check conservation (C Rousset) … … 108 110 CALL lim_thd_lac 109 111 CALL lim_var_glo2eqv ! only for info 110 111 IF(ln_ctl) THEN ! Control print112 113 IF(ln_ctl) THEN ! Control print 112 114 CALL prt_ctl_info(' ') 113 115 CALL prt_ctl_info(' - Cell values : ') … … 144 146 !- check conservation (C Rousset) 145 147 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 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 149 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 150 zchk_ft = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) - zchk_ft_b 148 151 149 zchk_v_i = ( glob_sum( SUM( v_i(:,:,:) , dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) * r1_rdtice152 zchk_v_i = ( glob_sum( SUM( v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b ) * r1_rdtice - zchk_fw 150 153 zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic ) 154 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 151 155 152 156 zchk_vmin = glob_min(v_i) … … 155 159 156 160 IF(lwp) THEN 157 IF ( ABS( zchk_v_i ) > 1.e- 5 ) WRITE(numout,*) 'violation volume [m3/day] (limitd_th) = ',(zchk_v_i * rday)161 IF ( ABS( zchk_v_i ) > 1.e-4 ) WRITE(numout,*) 'violation volume [kg/day] (limitd_th) = ',(zchk_v_i * rday) 158 162 IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limitd_th) = ',(zchk_smv * rday) 163 IF ( ABS( zchk_e_i ) > 1.e-2 ) WRITE(numout,*) 'violation enthalpy [1e9 J] (limitd_th) = ',(zchk_e_i) 159 164 IF ( zchk_vmin < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limitd_th) = ',(zchk_vmin * 1.e-3) 160 165 IF ( zchk_amax > amax+epsi10 ) WRITE(numout,*) 'violation a_i>amax (limitd_th) = ',zchk_amax … … 258 263 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 264 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)265 IF( a_i(ji,jj,jl) > epsi10 ) zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_o(ji,jj,jl) 261 266 END DO 262 267 END DO … … 302 307 ij = nind_j(ji) 303 308 ! 304 IF ( ( zht_i_o(ii,ij,jl) .GT. epsi10 ) .AND. &305 ( zht_i_o(ii,ij,jl+1) .GT. epsi10 )) THEN309 zhbnew(ii,ij,jl) = hi_max(jl) 310 IF ( old_a_i(ii,ij,jl) > epsi10 .AND. old_a_i(ii,ij,jl+1) > epsi10 ) THEN 306 311 !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 312 zslope = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / ( zht_i_o(ii,ij,jl+1) - zht_i_o(ii,ij,jl) ) 313 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + zslope * ( hi_max(jl) - zht_i_o(ii,ij,jl) ) 314 ELSEIF ( old_a_i(ii,ij,jl) > epsi10) THEN 312 315 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) 313 ELSEIF ( zht_i_o(ii,ij,jl+1).gt.epsi10) THEN316 ELSEIF ( old_a_i(ii,ij,jl+1) > epsi10) THEN 314 317 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1) 315 ELSE316 zhbnew(ii,ij,jl) = hi_max(jl)317 318 ENDIF 318 319 END DO … … 320 321 !- 4.2 Check that each zhbnew lies between adjacent values of ice thickness 321 322 DO ji = 1, nbrem 322 ! jl, ji323 323 ii = nind_i(ji) 324 324 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 325 IF( a_i(ii,ij,jl) > epsi10 .AND. ht_i(ii,ij,jl) >= zhbnew(ii,ij,jl) ) THEN 329 326 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 327 ELSEIF( a_i(ii,ij,jl+1) > epsi10 .AND. ht_i(ii,ij,jl+1) <= zhbnew(ii,ij,jl) ) THEN 333 328 zremap_flag(ii,ij) = 0 334 329 ENDIF 335 330 336 331 !- 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 332 IF( zhbnew(ii,ij,jl) > hi_max(jl+1) ) zremap_flag(ii,ij) = 0 333 IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0 334 END DO 335 348 336 END DO !jl 349 337 … … 354 342 DO jj = 1, jpj 355 343 DO ji = 1, jpi 356 IF 344 IF( zremap_flag(ji,jj) == 1 ) THEN 357 345 nbrem = nbrem + 1 358 346 nind_i(nbrem) = ji 359 347 nind_j(nbrem) = jj 360 348 ENDIF 361 END DO !ji362 END DO !jj349 END DO 350 END DO 363 351 364 352 !----------------------------------------------------------------------------------------------- … … 380 368 ENDIF 381 369 382 IF( zhbnew(ji,jj,kubnd) < hi_max(kubnd-1) ) 370 IF( zhbnew(ji,jj,kubnd) < hi_max(kubnd-1) ) zhbnew(ji,jj,kubnd) = hi_max(kubnd-1) 383 371 384 372 END DO !jj … … 444 432 DO jl = klbnd, kubnd 445 433 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) 434 g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), zremap_flag) 448 435 END DO 449 436 … … 493 480 nd = zdonor(ii,ij,jl) 494 481 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) 482 zdvice(ii,ij,jl) = g1(ii,ij,nd)*zx3 + g0(ii,ij,nd)*zx2 + zdaice(ii,ij,jl)*hL(ii,ij,nd) 497 483 498 484 END DO ! ji … … 511 497 ii = nind_i(ji) 512 498 ij = nind_j(ji) 513 IF ( ( a_i(ii,ij,1) > epsi10 ) .AND. ( ht_i(ii,ij,1) < hiclim )) THEN499 IF ( a_i(ii,ij,1) > epsi10 .AND. ht_i(ii,ij,1) < hiclim ) THEN 514 500 a_i(ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) / hiclim 515 501 ht_i(ii,ij,1) = hiclim 516 v_i(ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) !clem-useless517 502 ENDIF 518 503 END DO !ji … … 799 784 !-------------- 800 785 801 zdvsnow 786 zdvsnow = v_s(ii,ij,jl1) * zworka(ii,ij) 802 787 v_s(ii,ij,jl1) = v_s(ii,ij,jl1) - zdvsnow 803 788 v_s(ii,ij,jl2) = v_s(ii,ij,jl2) + zdvsnow … … 807 792 !-------------------- 808 793 809 zdesnow 794 zdesnow = e_s(ii,ij,1,jl1) * zworka(ii,ij) 810 795 e_s(ii,ij,1,jl1) = e_s(ii,ij,1,jl1) - zdesnow 811 796 e_s(ii,ij,1,jl2) = e_s(ii,ij,1,jl2) + zdesnow … … 815 800 !-------------- 816 801 817 zdo_aice 802 zdo_aice = oa_i(ii,ij,jl1) * zdaice(ii,ij,jl) 818 803 oa_i(ii,ij,jl1) = oa_i(ii,ij,jl1) - zdo_aice 819 804 oa_i(ii,ij,jl2) = oa_i(ii,ij,jl2) + zdo_aice … … 823 808 !-------------- 824 809 825 zdsm_vice 810 zdsm_vice = smv_i(ii,ij,jl1) * zworka(ii,ij) 826 811 smv_i(ii,ij,jl1) = smv_i(ii,ij,jl1) - zdsm_vice 827 812 smv_i(ii,ij,jl2) = smv_i(ii,ij,jl2) + zdsm_vice … … 831 816 !--------------------- 832 817 833 zdaTsf 818 zdaTsf = t_su(ii,ij,jl1) * zdaice(ii,ij,jl) 834 819 zaTsfn(ii,ij,jl1) = zaTsfn(ii,ij,jl1) - zdaTsf 835 820 zaTsfn(ii,ij,jl2) = zaTsfn(ii,ij,jl2) + zdaTsf … … 910 895 REAL(wp), POINTER, DIMENSION(:,:) :: vt_s_init, vt_s_final ! snow volume summed over categories 911 896 !!------------------------------------------------------------------ 897 !! clem 2014/04: be carefull, rebining does not conserve salt => the difference is taken into account in limupdate 912 898 913 899 CALL wrk_alloc( jpi,jpj,jpl, zdonor ) ! interger … … 1015 1001 1016 1002 !clem-change 1003 DO jj = 1, jpj 1004 DO ji = 1, jpi 1005 IF( a_i(ji,jj,jl+1) > epsi10 .AND. ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN 1006 ! 1007 zshiftflag = 1 1008 zdonor(ji,jj,jl) = jl + 1 1009 zdaice(ji,jj,jl) = a_i(ji,jj,jl+1) 1010 zdvice(ji,jj,jl) = v_i(ji,jj,jl+1) 1011 ENDIF 1012 END DO ! ji 1013 END DO ! jj 1014 1015 IF(lk_mpp) CALL mpp_max( zshiftflag ) 1016 1017 IF( zshiftflag == 1 ) THEN ! Shift ice between categories 1018 CALL lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice ) 1019 ! Reset shift parameters 1020 zdonor(:,:,jl) = 0 1021 zdaice(:,:,jl) = 0._wp 1022 zdvice(:,:,jl) = 0._wp 1023 ENDIF 1024 !clem-change 1025 1026 ! ! clem-change begin: why not doing that? 1017 1027 ! DO jj = 1, jpj 1018 1028 ! 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) 1029 ! IF( a_i(ji,jj,jl+1) > epsi10 .AND. ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN 1030 ! ht_i(ji,jj,jl+1) = hi_max(jl) + epsi10 1031 ! a_i (ji,jj,jl+1) = v_i(ji,jj,jl+1) / ht_i(ji,jj,jl+1) 1026 1032 ! ENDIF 1027 1033 ! END DO ! ji 1028 1034 ! 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 1035 ! clem-change end 1052 1036
Note: See TracChangeset
for help on using the changeset viewer.
