- Timestamp:
- 2015-03-24T18:35:00+01:00 (9 years ago)
- Location:
- trunk/NEMOGCM/NEMO
- Files:
-
- 18 edited
-
LIM_SRC_3/ice.F90 (modified) (4 diffs)
-
LIM_SRC_3/limcons.F90 (modified) (7 diffs)
-
LIM_SRC_3/limctl.F90 (modified) (1 diff)
-
LIM_SRC_3/limdiahsb.F90 (modified) (2 diffs)
-
LIM_SRC_3/limitd_me.F90 (modified) (15 diffs)
-
LIM_SRC_3/limsbc.F90 (modified) (2 diffs)
-
LIM_SRC_3/limthd.F90 (modified) (13 diffs)
-
LIM_SRC_3/limthd_dh.F90 (modified) (27 diffs)
-
LIM_SRC_3/limthd_dif.F90 (modified) (3 diffs)
-
LIM_SRC_3/limthd_lac.F90 (modified) (2 diffs)
-
LIM_SRC_3/limtrp.F90 (modified) (8 diffs)
-
LIM_SRC_3/limupdate1.F90 (modified) (4 diffs)
-
LIM_SRC_3/limupdate2.F90 (modified) (6 diffs)
-
LIM_SRC_3/limvar.F90 (modified) (11 diffs)
-
LIM_SRC_3/limwri.F90 (modified) (1 diff)
-
LIM_SRC_3/thd_ice.F90 (modified) (3 diffs)
-
OPA_SRC/BDY/bdyice_lim.F90 (modified) (3 diffs)
-
OPA_SRC/SBC/sbcice_lim.F90 (modified) (9 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
r5146 r5167 394 394 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_trp_smv !: transport of salt content 395 395 ! 396 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_heat_dhc !: snw/ice heat content variation [W/m2] 396 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_heat !: snw/ice heat content variation [W/m2] 397 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_smvi !: ice salt content variation [] 398 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_vice !: ice volume variation [m/s] 399 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_vsnw !: snw volume variation [m/s] 397 400 ! 398 401 !!---------------------------------------------------------------------- … … 455 458 ALLOCATE( t_s(jpi,jpj,nlay_s,jpl) , e_s(jpi,jpj,nlay_s,jpl) , STAT=ierr(ii) ) 456 459 ii = ii + 1 457 ALLOCATE( t_i(jpi,jpj,nlay_i +1,jpl) , e_i(jpi,jpj,nlay_i+1,jpl) , s_i(jpi,jpj,nlay_i+1,jpl) , STAT=ierr(ii) )460 ALLOCATE( t_i(jpi,jpj,nlay_i,jpl) , e_i(jpi,jpj,nlay_i,jpl) , s_i(jpi,jpj,nlay_i,jpl) , STAT=ierr(ii) ) 458 461 459 462 ! * Moments for advection … … 471 474 & STAT=ierr(ii) ) 472 475 ii = ii + 1 473 ALLOCATE( sxe (jpi,jpj,nlay_i +1,jpl) , sye (jpi,jpj,nlay_i+1,jpl) , sxxe(jpi,jpj,nlay_i+1,jpl) , &474 & syye(jpi,jpj,nlay_i +1,jpl) , sxye(jpi,jpj,nlay_i+1,jpl), STAT=ierr(ii) )476 ALLOCATE( sxe (jpi,jpj,nlay_i,jpl) , sye (jpi,jpj,nlay_i,jpl) , sxxe(jpi,jpj,nlay_i,jpl) , & 477 & syye(jpi,jpj,nlay_i,jpl) , sxye(jpi,jpj,nlay_i,jpl) , STAT=ierr(ii) ) 475 478 476 479 ! * Old values of global variables 477 480 ii = ii + 1 478 481 ALLOCATE( v_s_b (jpi,jpj,jpl) , v_i_b (jpi,jpj,jpl) , e_s_b(jpi,jpj,nlay_s,jpl) , & 479 & a_i_b (jpi,jpj,jpl) , smv_i_b(jpi,jpj,jpl) , e_i_b(jpi,jpj,nlay_i +1 ,jpl) ,&480 & oa_i_b (jpi,jpj,jpl) , u_ice_b(jpi,jpj) , v_ice_b(jpi,jpj) 482 & a_i_b (jpi,jpj,jpl) , smv_i_b(jpi,jpj,jpl) , e_i_b(jpi,jpj,nlay_i,jpl) , & 483 & oa_i_b (jpi,jpj,jpl) , u_ice_b(jpi,jpj) , v_ice_b(jpi,jpj) , STAT=ierr(ii) ) 481 484 482 485 ! * Ice thickness distribution variables … … 486 489 ! * Ice diagnostics 487 490 ii = ii + 1 488 ALLOCATE( diag_trp_vi(jpi,jpj), diag_trp_vs (jpi,jpj), diag_trp_ei (jpi,jpj), & 489 & diag_trp_es(jpi,jpj), diag_trp_smv(jpi,jpj), diag_heat_dhc(jpi,jpj), STAT=ierr(ii) ) 491 ALLOCATE( diag_trp_vi(jpi,jpj), diag_trp_vs (jpi,jpj), diag_trp_ei(jpi,jpj), & 492 & diag_trp_es(jpi,jpj), diag_trp_smv(jpi,jpj), diag_heat (jpi,jpj), & 493 & diag_smvi (jpi,jpj), diag_vice (jpi,jpj), diag_vsnw (jpi,jpj), STAT=ierr(ii) ) 490 494 491 495 ice_alloc = MAXVAL( ierr(:) ) -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90
r5123 r5167 22 22 USE lib_mpp ! MPP library 23 23 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 24 USE sbc_oce , ONLY : sfx ! Surface boundary condition: ocean fields 24 25 25 26 IMPLICIT NONE … … 30 31 PUBLIC lim_cons_check 31 32 PUBLIC lim_cons_hsm 33 PUBLIC lim_cons_final 32 34 33 35 !!---------------------------------------------------------------------- … … 167 169 REAL(wp) :: zvi, zsmv, zei, zfs, zfw, zft 168 170 REAL(wp) :: zvmin, zamin, zamax 169 REAL(wp) :: zconv 170 171 zconv = 1.e-9 171 REAL(wp) :: zvtrp, zetrp 172 REAL(wp), PARAMETER :: zconv = 1.e-9 172 173 173 174 IF( icount == 0 ) THEN … … 187 188 zvi_b = glob_sum( SUM( v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * e12t(:,:) * tmask(:,:,1) ) 188 189 189 zsmv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) )190 zsmv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) * rhoic ) 190 191 191 192 zei_b = glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) + & … … 210 211 & * e12t(:,:) * tmask(:,:,1) ) - zvi_b ) * r1_rdtice - zfw 211 212 212 zsmv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) ) - zsmv_b ) * r1_rdtice + ( zfs * r1_rhoic )213 zsmv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) * rhoic ) - zsmv_b ) * r1_rdtice + zfs 213 214 214 215 zei = glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) + & … … 216 217 & ) * e12t(:,:) * tmask(:,:,1) * zconv ) * r1_rdtice - zei_b * r1_rdtice + zft 217 218 219 zvtrp = glob_sum( ( diag_trp_vi * rhoic + diag_trp_vs * rhosn ) * e12t(:,:) * tmask(:,:,1) ) 220 zetrp = glob_sum( ( diag_trp_ei + diag_trp_es ) * e12t(:,:) * tmask(:,:,1) * zconv ) 218 221 zvmin = glob_min( v_i ) 219 222 zamax = glob_max( SUM( a_i, dim=3 ) ) 220 223 zamin = glob_min( a_i ) 224 221 225 222 226 IF(lwp) THEN 223 IF ( ABS( zvi ) > 1.e-4) WRITE(numout,*) 'violation volume [kg/day] (',cd_routine,') = ',(zvi * rday)224 IF ( ABS( zsmv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (',cd_routine,') = ',(zsmv * rday)225 IF ( ABS( zei ) > 1.e-4) WRITE(numout,*) 'violation enthalpy [GW] (',cd_routine,') = ',(zei)226 IF ( zvmin < -epsi10 ) WRITE(numout,*) 'violation v_i<0 [m] (',cd_routine,') = ',(zvmin)227 IF ( ABS( zvi * rday ) > 0.5 * 1.e9 ) WRITE(numout,*) 'violation volume [kg/day] (',cd_routine,') = ',(zvi * rday) 228 IF ( ABS( zsmv * rday ) > 5. * 1.e9 ) WRITE(numout,*) 'violation saline [psu*kg/day] (',cd_routine,') = ',(zsmv * rday) 229 IF ( ABS( zei ) > 2. * 1.e9 ) WRITE(numout,*) 'violation enthalpy [GW] (',cd_routine,') = ',(zei) 230 IF ( zvmin < -epsi10 ) WRITE(numout,*) 'violation v_i<0 [m] (',cd_routine,') = ',(zvmin) 227 231 IF( cd_routine /= 'limtrp' .AND. cd_routine /= 'limitd_me' .AND. zamax > rn_amax+epsi10 ) THEN 228 WRITE(numout,*) 'violation a_i>amax (',cd_routine,') = ',zamax232 WRITE(numout,*) 'violation a_i>amax (',cd_routine,') = ',zamax 229 233 ENDIF 230 IF ( zamin < -epsi10 ) WRITE(numout,*) 'violation a_i<0 (',cd_routine,') = ',zamin 234 IF ( zamin < -epsi10 ) WRITE(numout,*) 'violation a_i<0 (',cd_routine,') = ',zamin 235 IF( cd_routine == 'limtrp' .AND. ABS( zvtrp * rday ) > 0.5*1.e9 ) THEN 236 WRITE(numout,*) 'violation vtrp [kg/day] (',cd_routine,') = ',(zvtrp * rday) 237 WRITE(numout,*) 'violation etrp [GW] (',cd_routine,') = ',(zetrp ) 238 ENDIF 231 239 ENDIF 232 240 … … 234 242 235 243 END SUBROUTINE lim_cons_hsm 244 245 SUBROUTINE lim_cons_final( cd_routine ) 246 CHARACTER(len=*), INTENT(in) :: cd_routine ! name of the routine 247 REAL(wp) :: zhfx, zsfx, zvfx 248 REAL(wp), PARAMETER :: zconv = 1.e-9 249 250 zhfx = glob_sum( ( hfx_in - hfx_out - diag_heat - diag_trp_ei - diag_trp_es - hfx_sub ) * e12t(:,:) * tmask(:,:,1) * zconv ) 251 zsfx = glob_sum( ( sfx + diag_smvi ) * e12t(:,:) * tmask(:,:,1) ) * rday 252 zvfx = glob_sum( ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e12t(:,:) * tmask(:,:,1) ) * rday 253 254 ! if error > 1 mm / 100 years over the Arctic Basin 255 IF( ABS( zvfx ) > 0.5 * 1.e9 ) WRITE(numout,*) 'violation vfx [kg/day] (',cd_routine,') = ',(zvfx) 256 ! if error > 1 mm / 100 years over the Arctic Basin (ice with latent heat = 3e6 J/kg) 257 IF( ABS( zhfx ) > 2. * 1.e9 ) WRITE(numout,*) 'violation hfx [GW] (',cd_routine,') = ',(zhfx) 258 ! if error > 1 mm / 100 years over the Arctic Basin (ice of salinity = 10 pss) 259 IF( ABS( zsfx ) > 5. * 1.e9 ) WRITE(numout,*) 'violation sfx [psu*kg/day] (',cd_routine,') = ',(zsfx) 260 261 END SUBROUTINE lim_cons_final 236 262 237 263 #else -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limctl.F90
r5125 r5167 419 419 WRITE(numout,*) ' hfx_in : ', hfx_in(ji,jj) 420 420 WRITE(numout,*) ' hfx_out : ', hfx_out(ji,jj) 421 WRITE(numout,*) ' dhc : ', diag_heat _dhc(ji,jj)421 WRITE(numout,*) ' dhc : ', diag_heat(ji,jj) 422 422 WRITE(numout,*) 423 423 WRITE(numout,*) ' hfx_dyn : ', hfx_dyn(ji,jj) -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90
r5123 r5167 115 115 zbg_ihc = glob_sum( et_i(:,:) * e12t(:,:) * 1.e-20 ) ! ice heat content [1.e20 J] 116 116 zbg_shc = glob_sum( et_s(:,:) * e12t(:,:) * 1.e-20 ) ! snow heat content [1.e20 J] 117 zbg_hfx_dhc = glob_sum( diag_heat _dhc(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W]117 zbg_hfx_dhc = glob_sum( diag_heat(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 118 118 zbg_hfx_spr = glob_sum( hfx_spr(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 119 119 … … 245 245 WRITE(numout,*) '~~~~~~~~~~~~' 246 246 ENDIF 247 248 ! ---------------------------------- !249 ! 2 - initial conservation variables !250 ! ---------------------------------- !251 !frc_vol = 0._wp ! volume trend due to forcing252 !frc_sal = 0._wp ! salt content - - - -253 !bg_grme = 0._wp ! ice growth + melt volume trend254 247 ! 255 248 CALL lim_diahsb_rst( nstart, 'READ' ) !* read or initialize all required files -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
r5134 r5167 154 154 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 155 155 156 CALL lim_var_glo2eqv ! equivalent variables, requested for rafting 156 157 !-----------------------------------------------------------------------------! 157 158 ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons … … 235 236 ! Reduce the closing rate if more than 100% of the open water 236 237 ! would be removed. Reduce the opening rate proportionately. 237 IF ( ato_i(ji,jj) > epsi10 .AND. athorn(ji,jj,0) > 0.0 ) THEN 238 za = athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice 239 IF ( za > ato_i(ji,jj)) THEN 240 zfac = ato_i(ji,jj) / za 241 closing_gross(ji,jj) = closing_gross(ji,jj) * zfac 242 opning(ji,jj) = opning(ji,jj) * zfac 243 ENDIF 238 za = athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice 239 IF( za > epsi20 ) THEN 240 zfac = MIN( 1._wp, ato_i(ji,jj) / za ) 241 closing_gross(ji,jj) = closing_gross(ji,jj) * zfac 242 opning (ji,jj) = opning (ji,jj) * zfac 244 243 ENDIF 245 244 … … 251 250 ! Reduce the closing rate if more than 100% of any ice category 252 251 ! would be removed. Reduce the opening rate proportionately. 253 254 252 DO jl = 1, jpl 255 253 DO jj = 1, jpj 256 254 DO ji = 1, jpi 257 IF ( a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0._wp )THEN 258 za = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice 259 IF ( za > a_i(ji,jj,jl) ) THEN 260 zfac = a_i(ji,jj,jl) / za 261 closing_gross(ji,jj) = closing_gross(ji,jj) * zfac 262 opning (ji,jj) = opning (ji,jj) * zfac 263 ENDIF 255 za = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice 256 IF( za > epsi20 ) THEN 257 zfac = MIN( 1._wp, a_i(ji,jj,jl) / za ) 258 closing_gross(ji,jj) = closing_gross(ji,jj) * zfac 259 opning (ji,jj) = opning (ji,jj) * zfac 264 260 ENDIF 265 261 END DO … … 369 365 370 366 ! updates 371 CALL lim_var_glo2eqv372 CALL lim_var_zapsmall373 367 CALL lim_var_agg( 1 ) 374 368 … … 377 371 !-----------------------------------------------------------------------------! 378 372 IF(ln_ctl) THEN 373 CALL lim_var_glo2eqv 374 379 375 CALL prt_ctl_info(' ') 380 376 CALL prt_ctl_info(' - Cell values : ') … … 531 527 DO jj = 2, jpjm1 532 528 DO ji = 2, jpim1 533 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > epsi10) THEN ! ice is present529 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp) THEN 534 530 zworka(ji,jj) = ( 4.0 * strength(ji,jj) & 535 531 & + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) & … … 566 562 DO jj = 1, jpj - 1 567 563 DO ji = 1, jpi - 1 568 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > epsi10) THEN ! ice is present564 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp) THEN 569 565 numts_rm = 1 ! number of time steps for the running mean 570 566 IF ( strp1(ji,jj) > 0.0 ) numts_rm = numts_rm + 1 … … 637 633 638 634 Gsum(:,:,-1) = 0._wp 639 640 DO jj = 1, jpj 641 DO ji = 1, jpi 642 IF( ato_i(ji,jj) > epsi10 ) THEN ; Gsum(ji,jj,0) = ato_i(ji,jj) 643 ELSE ; Gsum(ji,jj,0) = 0._wp 644 ENDIF 645 END DO 646 END DO 635 Gsum(:,:,0 ) = ato_i(:,:) 647 636 648 637 ! for each value of h, you have to add ice concentration then 649 638 DO jl = 1, jpl 650 DO jj = 1, jpj 651 DO ji = 1, jpi 652 IF( a_i(ji,jj,jl) > epsi10 ) THEN ; Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) + a_i(ji,jj,jl) 653 ELSE ; Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) 654 ENDIF 655 END DO 656 END DO 639 Gsum(:,:,jl) = Gsum(:,:,jl-1) + a_i(:,:,jl) 657 640 END DO 658 641 … … 828 811 LOGICAL, PARAMETER :: l_conservation_check = .true. ! if true, check conservation (useful for debugging) 829 812 ! 830 LOGICAL :: neg_ato_i ! flag for ato_i(i,j) < -puny831 LOGICAL :: large_afrac ! flag for afrac > 1832 LOGICAL :: large_afrft ! flag for afrac > 1833 813 INTEGER :: ji, jj, jl, jl1, jl2, jk ! dummy loop indices 834 814 INTEGER :: ij ! horizontal index, combines i and j loops … … 898 878 ! 1) Compute change in open water area due to closing and opening. 899 879 !------------------------------------------------------------------------------- 900 901 neg_ato_i = .false.902 903 880 DO jj = 1, jpj 904 881 DO ji = 1, jpi 905 882 ato_i(ji,jj) = ato_i(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice & 906 883 & + opning(ji,jj) * rdt_ice 907 IF ( ato_i(ji,jj) < -epsi10 ) THEN908 neg_ato_i = .TRUE.909 ELSEIF( ato_i(ji,jj) < 0._wp ) THEN ! roundoff error884 IF ( ato_i(ji,jj) < -epsi10 ) THEN ! there is a bug 885 IF(lwp) WRITE(numout,*) 'Ridging error: ato_i < 0 -- ato_i : ',ato_i(ji,jj) 886 ELSEIF( ato_i(ji,jj) < 0._wp ) THEN ! roundoff error 910 887 ato_i(ji,jj) = 0._wp 911 888 ENDIF 912 889 END DO 913 890 END DO 914 915 ! if negative open water area alert it916 IF( neg_ato_i .AND. lwp ) THEN ! there is a bug917 DO jj = 1, jpj918 DO ji = 1, jpi919 IF( ato_i(ji,jj) < -epsi10 ) THEN920 WRITE(numout,*) ''921 WRITE(numout,*) 'Ridging error: ato_i < 0'922 WRITE(numout,*) 'ato_i : ', ato_i(ji,jj)923 ENDIF924 END DO925 END DO926 ENDIF927 891 928 892 !----------------------------------------------------------------- 929 893 ! 2) Save initial state variables 930 894 !----------------------------------------------------------------- 931 932 DO jl = 1, jpl 933 aicen_init(:,:,jl) = a_i(:,:,jl) 934 vicen_init(:,:,jl) = v_i(:,:,jl) 935 vsnwn_init(:,:,jl) = v_s(:,:,jl) 936 ! 937 smv_i_init(:,:,jl) = smv_i(:,:,jl) 938 oa_i_init (:,:,jl) = oa_i (:,:,jl) 939 END DO 940 941 esnwn_init(:,:,:) = e_s(:,:,1,:) 942 943 DO jl = 1, jpl 944 DO jk = 1, nlay_i 945 eicen_init(:,:,jk,jl) = e_i(:,:,jk,jl) 946 END DO 947 END DO 895 aicen_init(:,:,:) = a_i (:,:,:) 896 vicen_init(:,:,:) = v_i (:,:,:) 897 vsnwn_init(:,:,:) = v_s (:,:,:) 898 smv_i_init(:,:,:) = smv_i(:,:,:) 899 oa_i_init (:,:,:) = oa_i (:,:,:) 900 esnwn_init(:,:,:) = e_s (:,:,1,:) 901 eicen_init(:,:,:,:) = e_i (:,:,:,:) 948 902 949 903 ! … … 972 926 END DO 973 927 974 large_afrac = .false.975 large_afrft = .false.976 977 928 DO ij = 1, icells 978 929 ji = indxi(ij) … … 1000 951 afrft(ji,jj) = arft1(ji,jj) / aicen_init(ji,jj,jl1) !rafting 1001 952 1002 IF (afrac(ji,jj) > kamax + epsi10) THEN !riging1003 large_afrac = .true.1004 ELSEIF (afrac(ji,jj) > kamax) THEN! roundoff error953 IF( afrac(ji,jj) > kamax + epsi10 ) THEN ! there is a bug 954 IF(lwp) WRITE(numout,*) ' ardg > a_i -- ardg, aicen_init : ', ardg1(ji,jj), aicen_init(ji,jj,jl1) 955 ELSEIF( afrac(ji,jj) > kamax ) THEN ! roundoff error 1005 956 afrac(ji,jj) = kamax 1006 957 ENDIF 1007 IF (afrft(ji,jj) > kamax + epsi10) THEN !rafting 1008 large_afrft = .true. 1009 ELSEIF (afrft(ji,jj) > kamax) THEN ! roundoff error 958 959 IF( afrft(ji,jj) > kamax + epsi10 ) THEN ! there is a bug 960 IF(lwp) WRITE(numout,*) ' arft > a_i -- arft, aicen_init : ', arft1(ji,jj), aicen_init(ji,jj,jl1) 961 ELSEIF( afrft(ji,jj) > kamax) THEN ! roundoff error 1010 962 afrft(ji,jj) = kamax 1011 963 ENDIF … … 1022 974 esrdg(ji,jj) = esnwn_init(ji,jj,jl1) * afrac(ji,jj) 1023 975 srdg1(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) 1024 srdg2(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) !! MV HC 2014 this line seems useless1025 976 1026 977 ! rafting volumes, heat contents ... … … 1050 1001 1051 1002 sfx_dyn(ji,jj) = sfx_dyn(ji,jj) - smsw(ji,jj) * rhoic * r1_rdtice 1052 wfx_dyn(ji,jj) = wfx_dyn(ji,jj) - vsw (ji,jj) * rhoic * r1_rdtice ! gurvan:increase in ice volume du to seawater frozen in voids1003 wfx_dyn(ji,jj) = wfx_dyn(ji,jj) - vsw (ji,jj) * rhoic * r1_rdtice ! increase in ice volume du to seawater frozen in voids 1053 1004 1054 1005 !------------------------------------ … … 1134 1085 ENDIF 1135 1086 1136 IF( large_afrac .AND. lwp ) THEN ! there is a bug1137 DO ij = 1, icells1138 ji = indxi(ij)1139 jj = indxj(ij)1140 IF( afrac(ji,jj) > kamax + epsi10 ) THEN1141 WRITE(numout,*) ''1142 WRITE(numout,*) ' ardg > a_i'1143 WRITE(numout,*) ' ardg, aicen_init : ', ardg1(ji,jj), aicen_init(ji,jj,jl1)1144 ENDIF1145 END DO1146 ENDIF1147 IF( large_afrft .AND. lwp ) THEN ! there is a bug1148 DO ij = 1, icells1149 ji = indxi(ij)1150 jj = indxj(ij)1151 IF( afrft(ji,jj) > kamax + epsi10 ) THEN1152 WRITE(numout,*) ''1153 WRITE(numout,*) ' arft > a_i'1154 WRITE(numout,*) ' arft, aicen_init : ', arft1(ji,jj), aicen_init(ji,jj,jl1)1155 ENDIF1156 END DO1157 ENDIF1158 1159 1087 !------------------------------------------------------------------------------- 1160 1088 ! 4) Add area, volume, and energy of new ridge to each category jl2 -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
r5146 r5167 42 42 USE domvvl ! Variable volume 43 43 USE limctl 44 USE limcons 44 45 45 46 IMPLICIT NONE … … 227 228 ENDIF 228 229 229 IF( ln_icectl ) CALL lim_prt( kt, iiceprt, jiceprt, 3, ' - Final state lim_sbc - ' ) ! control print 230 ! conservation test 231 IF( ln_limdiahsb ) CALL lim_cons_final( 'limsbc' ) 232 233 ! control prints 234 IF( ln_icectl ) CALL lim_prt( kt, iiceprt, jiceprt, 3, ' - Final state lim_sbc - ' ) 230 235 231 236 IF(ln_ctl) THEN -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
r5146 r5167 94 94 REAL(wp), POINTER, DIMENSION(:,:) :: zqsr, zqns 95 95 !!------------------------------------------------------------------- 96 CALL wrk_alloc( jpi, 96 CALL wrk_alloc( jpi,jpj, zqsr, zqns ) 97 97 98 98 IF( nn_timing == 1 ) CALL timing_start('limthd') … … 101 101 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 102 102 103 CALL lim_var_glo2eqv 103 104 !------------------------------------------------------------------------! 104 105 ! 1) Initialization of some variables ! … … 209 210 ! Net heat flux on top of ice-ocean [W.m-2] 210 211 ! ----------------------------------------- 211 ! First step here :heat flux at the ocean surface + precip212 ! Second step below : heat flux at the ice surface (after limthd_dif)212 ! heat flux at the ocean surface + precip 213 ! + heat flux at the ice surface 213 214 hfx_in(ji,jj) = hfx_in(ji,jj) & 214 215 ! heat flux above the ocean … … 216 217 ! latent heat of precip (note that precip is included in qns but not in qns_ice) 217 218 & + ( 1._wp - pfrld(ji,jj) ) * sprecip(ji,jj) * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rt0 ) - lfus ) & 218 & + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rt0 ) 219 & + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rt0 ) & 220 ! heat flux above the ice 221 & + SUM( a_i_b(ji,jj,:) * ( qns_ice(ji,jj,:) + qsr_ice(ji,jj,:) ) ) 219 222 220 223 ! ----------------------------------------------------------------------------- … … 226 229 hfx_out(ji,jj) = hfx_out(ji,jj) & 227 230 ! Non solar heat flux received by the ocean 228 & + pfrld(ji,jj) * qns(ji,jj) &231 & + pfrld(ji,jj) * zqns(ji,jj) & 229 232 ! latent heat of precip (note that precip is included in qns but not in qns_ice) 230 233 & + ( pfrld(ji,jj)**rn_betas - pfrld(ji,jj) ) * sprecip(ji,jj) & … … 311 314 ! --- lateral melting if monocat --- ! 312 315 !------------------------------------! 313 IF ( ( ( nn_monocat == 1 ) .OR. ( nn_monocat == 4 ) ) .AND. ( jpl == 1 )) THEN316 IF ( ( nn_monocat == 1 .OR. nn_monocat == 4 ) .AND. jpl == 1 ) THEN 314 317 CALL lim_thd_lam( 1, nbpb ) 315 318 END IF … … 324 327 ENDIF 325 328 ! 326 END DO 329 END DO !jl 327 330 328 331 !------------------------------------------------------------------------------! … … 358 361 ! Change thickness to volume 359 362 !---------------------------------- 360 CALL lim_var_eqv2glo 363 v_i(:,:,:) = ht_i(:,:,:) * a_i(:,:,:) 364 v_s(:,:,:) = ht_s(:,:,:) * a_i(:,:,:) 365 smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:) 361 366 362 367 CALL lim_var_zapsmall … … 399 404 ! 400 405 ! 401 CALL wrk_dealloc( jpi, jpj, zqsr, zqns )402 403 406 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 407 408 CALL wrk_dealloc( jpi,jpj, zqsr, zqns ) 409 404 410 !------------------------------------------------------------------------------| 405 411 ! 6) Transport of ice between thickness categories. | 406 412 !------------------------------------------------------------------------------| 413 ! Given thermodynamic growth rates, transport ice between thickness categories. 407 414 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 408 415 409 ! Given thermodynamic growth rates, transport ice between thickness categories. 410 IF( jpl > 1 ) CALL lim_itd_th_rem( 1, jpl, kt ) 411 ! 412 CALL lim_var_glo2eqv ! only for info 413 CALL lim_var_agg(1) 416 IF( jpl > 1 ) CALL lim_itd_th_rem( 1, jpl, kt ) 414 417 415 418 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 419 416 420 !------------------------------------------------------------------------------| 417 421 ! 7) Add frazil ice growing in leads. 418 422 !------------------------------------------------------------------------------| 419 423 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 424 420 425 CALL lim_thd_lac 421 CALL lim_var_glo2eqv ! only for info422 426 423 ! conservation test424 427 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 425 428 426 IF(ln_ctl) THEN ! Control print 429 ! Control print 430 IF(ln_ctl) THEN 431 CALL lim_var_glo2eqv 432 427 433 CALL prt_ctl_info(' ') 428 434 CALL prt_ctl_info(' - Cell values : ') … … 503 509 REAL(wp) :: zhi_bef ! ice thickness before thermo 504 510 REAL(wp) :: zdh_mel, zda_mel ! net melting 505 REAL(wp) :: zv ! ice volume511 REAL(wp) :: zvi, zvs ! ice/snow volumes 506 512 507 513 DO ji = kideb, kiut 508 514 zdh_mel = MIN( 0._wp, dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) ) 509 IF( zdh_mel < 0._wp ) THEN 510 zv = a_i_1d(ji) * ht_i_1d(ji) 515 IF( zdh_mel < 0._wp .AND. a_i_1d(ji) > 0._wp ) THEN 516 zvi = a_i_1d(ji) * ht_i_1d(ji) 517 zvs = a_i_1d(ji) * ht_s_1d(ji) 511 518 ! lateral melting = concentration change 512 519 zhi_bef = ht_i_1d(ji) - zdh_mel 513 zda_mel = a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi10 ) ) 514 a_i_1d(ji) = MAX( 0._wp, a_i_1d(ji) + zda_mel ) 515 ! adjust thickness 516 rswitch = MAX( 0._wp , SIGN( 1._wp , a_i_1d(ji) - epsi20 ) ) 517 ht_i_1d(ji) = rswitch * zv / MAX( a_i_1d(ji), epsi20 ) 520 rswitch = MAX( 0._wp , SIGN( 1._wp , zhi_bef - epsi20 ) ) 521 zda_mel = rswitch * a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi20 ) ) 522 a_i_1d(ji) = MAX( epsi20, a_i_1d(ji) + zda_mel ) 523 ! adjust thickness 524 ht_i_1d(ji) = zvi / a_i_1d(ji) 525 ht_s_1d(ji) = zvs / a_i_1d(ji) 518 526 ! retrieve total concentration 519 527 at_i_1d(ji) = a_i_1d(ji) … … 676 684 INTEGER :: ios ! Local integer output status for namelist read 677 685 NAMELIST/namicethd/ rn_hnewice, ln_frazil, rn_maxfrazb, rn_vfrazb, rn_Cfrazb, & 678 & rn_himin, parsub,rn_betas, rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon, &686 & rn_himin, rn_betas, rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon, & 679 687 & nn_monocat, ln_it_qnsice 680 688 !!------------------------------------------------------------------- … … 700 708 ENDIF 701 709 702 IF( lk_cpl .AND. parsub /= 0.0 ) CALL ctl_stop( 'In coupled mode, use parsub = 0. or send dqla' )703 710 ! 704 711 IF(lwp) THEN ! control print … … 712 719 WRITE(numout,*)' minimum ice thickness rn_himin = ', rn_himin 713 720 WRITE(numout,*)' numerical carac. of the scheme for diffusion in ice ' 714 WRITE(numout,*)' switch for snow sublimation (=1) or not (=0) parsub = ', parsub715 721 WRITE(numout,*)' coefficient for ice-lead partition of snowfall rn_betas = ', rn_betas 716 722 WRITE(numout,*)' extinction radiation parameter in sea ice rn_kappa_i = ', rn_kappa_i -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
r5134 r5167 86 86 REAL(wp) :: zsstK ! SST in Kelvin 87 87 88 REAL(wp), POINTER, DIMENSION(:) :: zh_s ! snow layer thickness89 88 REAL(wp), POINTER, DIMENSION(:) :: zqprec ! energy of fallen snow (J.m-3) 90 89 REAL(wp), POINTER, DIMENSION(:) :: zq_su ! heat for surface ablation (J.m-2) … … 118 117 END SELECT 119 118 120 CALL wrk_alloc( jpij, z h_s, zqprec, zq_su, zq_bo, zf_tt, zq_rema )119 CALL wrk_alloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema ) 121 120 CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s ) 122 121 CALL wrk_alloc( jpij, nlay_i+1, zdeltah, zh_i ) … … 129 128 zq_rema(:) = 0._wp 130 129 131 zh_s (:) = 0._wp132 130 zdh_s_pre(:) = 0._wp 133 131 zdh_s_mel(:) = 0._wp … … 140 138 icount (:) = 0 141 139 140 ! Initialize enthalpy at nlay_i+1 141 DO ji = kideb, kiut 142 q_i_1d(ji,nlay_i+1) = 0._wp 143 END DO 144 142 145 ! initialize layer thicknesses and enthalpies 143 146 h_i_old (:,0:nlay_i+1) = 0._wp … … 155 158 ! 156 159 DO ji = kideb, kiut 157 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) )158 ztmelts = rswitch * rt0 + ( 1._wp - rswitch ) * rt0159 160 160 zfdum = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 161 161 zf_tt(ji) = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji) 162 162 163 zq_su (ji) = MAX( 0._wp, zfdum * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - ztmelts) )163 zq_su (ji) = MAX( 0._wp, zfdum * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) ) 164 164 zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice ) 165 165 END DO … … 187 187 !------------------------------------------------------------! 188 188 ! 189 DO ji = kideb, kiut190 zh_s(ji) = ht_s_1d(ji) * r1_nlay_s191 END DO192 !193 189 DO jk = 1, nlay_s 194 190 DO ji = kideb, kiut 195 zqh_s(ji) = zqh_s(ji) + q_s_1d(ji,jk) * zh_s(ji)191 zqh_s(ji) = zqh_s(ji) + q_s_1d(ji,jk) * ht_s_1d(ji) * r1_nlay_s 196 192 END DO 197 193 END DO … … 222 218 ! Martin Vancoppenolle, December 2006 223 219 220 zdeltah(:,:) = 0._wp 224 221 DO ji = kideb, kiut 225 222 !----------- … … 236 233 ! mass flux, <0 237 234 wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice 238 ! update thickness239 ht_s_1d (ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_pre(ji) )240 235 241 236 !--------------------- … … 243 238 !--------------------- 244 239 ! thickness change 245 IF( zdh_s_pre(ji) > 0._wp ) THEN246 240 rswitch = MAX( 0._wp , SIGN( 1._wp , zqprec(ji) - epsi20 ) ) 247 zd h_s_mel (ji) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 )248 zd h_s_mel (ji) = MAX( - zdh_s_pre(ji), zdh_s_mel(ji) ) ! bound melting241 zdeltah (ji,1) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 ) 242 zdeltah (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) ) ! bound melting 249 243 ! heat used to melt snow (W.m-2, >0) 250 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zd h_s_mel(ji) * a_i_1d(ji) * zqprec(ji) * r1_rdtice244 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji) * r1_rdtice 251 245 ! snow melting only = water into the ocean (then without snow precip), >0 252 wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_mel(ji) * r1_rdtice 253 254 ! updates available heat + thickness 255 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdh_s_mel(ji) * zqprec(ji) ) 256 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_mel(ji) ) 257 zh_s (ji) = ht_s_1d(ji) * r1_nlay_s 258 259 ENDIF 260 END DO 261 262 ! If heat still available, then melt more snow 263 zdeltah(:,:) = 0._wp ! important 246 wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice 247 ! updates available heat + precipitations after melting 248 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,1) * zqprec(ji) ) 249 zdh_s_pre (ji) = zdh_s_pre(ji) + zdeltah(ji,1) 250 251 ! update thickness 252 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_pre(ji) ) 253 END DO 254 255 ! If heat still available (zq_su > 0), then melt more snow 256 zdeltah(:,:) = 0._wp 264 257 DO jk = 1, nlay_s 265 258 DO ji = kideb, kiut … … 268 261 rswitch = rswitch * ( MAX( 0._wp, SIGN( 1._wp, q_s_1d(ji,jk) - epsi20 ) ) ) 269 262 zdeltah (ji,jk) = - rswitch * zq_su(ji) / MAX( q_s_1d(ji,jk), epsi20 ) 270 zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji) ) ! bound melting263 zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - ht_s_1d(ji) ) ! bound melting 271 264 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk) 272 265 ! heat used to melt snow(W.m-2, >0) … … 274 267 ! snow melting only = water into the ocean (then without snow precip) 275 268 wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 276 277 269 ! updates available heat + thickness 278 270 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) ) 279 271 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdeltah(ji,jk) ) 280 281 272 END DO 282 273 END DO … … 286 277 !---------------------- 287 278 ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates 288 ! clem comment: not counted in mass exchange in limsbc since this is an exchange with atm. (not ocean)279 ! clem comment: not counted in mass/heat exchange in limsbc since this is an exchange with atm. (not ocean) 289 280 ! clem comment: ice should also sublimate 281 zdeltah(:,:) = 0._wp 290 282 IF( lk_cpl ) THEN 291 283 ! coupled mode: sublimation already included in emp_ice (to do in limsbc_ice) … … 294 286 ! forced mode: snow thickness change due to sublimation 295 287 DO ji = kideb, kiut 296 zdh_s_sub(ji) = MAX( - ht_s_1d(ji) , - parsub *qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice )288 zdh_s_sub(ji) = MAX( - ht_s_1d(ji) , - qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice ) 297 289 ! Heat flux by sublimation [W.m-2], < 0 298 290 ! sublimate first snow that had fallen, then pre-existing snow 299 zcoeff = ( MAX( zdh_s_sub(ji), - MAX( 0._wp, zdh_s_pre(ji) + zdh_s_mel(ji) ) ) * zqprec(ji) + & 300 & ( zdh_s_sub(ji) - MAX( zdh_s_sub(ji), - MAX( 0._wp, zdh_s_pre(ji) + zdh_s_mel(ji) ) ) ) * q_s_1d(ji,1) ) & 301 & * a_i_1d(ji) * r1_rdtice 302 hfx_sub_1d(ji) = hfx_sub_1d(ji) + zcoeff 291 zdeltah(ji,1) = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) ) 292 hfx_sub_1d(ji) = hfx_sub_1d(ji) + ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * q_s_1d(ji,1) & 293 & ) * a_i_1d(ji) * r1_rdtice 303 294 ! Mass flux by sublimation 304 295 wfx_sub_1d(ji) = wfx_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice 305 296 ! new snow thickness 306 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_sub(ji) ) 297 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_sub(ji) ) 298 ! update precipitations after sublimation and correct sublimation 299 zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1) 300 zdh_s_sub(ji) = zdh_s_sub(ji) - zdeltah(ji,1) 307 301 END DO 308 302 ENDIF … … 310 304 ! --- Update snow diags --- ! 311 305 DO ji = kideb, kiut 312 dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji) 313 zh_s(ji) = ht_s_1d(ji) * r1_nlay_s 314 END DO ! ji 306 dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji) 307 END DO 315 308 316 309 !------------------------------------------- … … 323 316 rswitch = MAX( 0._wp , SIGN( 1._wp, ht_s_1d(ji) - epsi20 ) ) 324 317 q_s_1d(ji,jk) = rswitch / MAX( ht_s_1d(ji), epsi20 ) * & 325 & ( ( MAX( 0._wp, dh_s_tot(ji) )) * zqprec(ji) + &326 & ( - MAX( 0._wp, dh_s_tot(ji) ) + ht_s_1d(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) )318 & ( ( zdh_s_pre(ji) ) * zqprec(ji) + & 319 & ( ht_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) ) 327 320 zq_s(ji) = zq_s(ji) + q_s_1d(ji,jk) 328 321 END DO … … 334 327 zdeltah(:,:) = 0._wp ! important 335 328 DO jk = 1, nlay_i 336 DO ji = kideb, kiut 337 zEi = - q_i_1d(ji,jk) * r1_rhoic 338 339 ztmelts = - tmut * s_i_1d(ji,jk) + rt0 329 DO ji = kideb, kiut 330 zEi = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of layer k [J/kg, <0] 331 332 ztmelts = - tmut * s_i_1d(ji,jk) + rt0 ! Melting point of layer k [K] 340 333 341 334 zEw = rcp * ( ztmelts - rt0 ) ! Specific enthalpy of resulting meltwater [J/kg, <0] … … 343 336 zdE = zEi - zEw ! Specific enthalpy difference < 0 344 337 345 zfmdt = - zq_su(ji) / zdE ! Mass flux to the ocean [kg/m2, >0] 338 IF( zdE < 0._wp ) THEN 339 zfmdt = - zq_su(ji) / zdE ! Mass flux to the ocean [kg/m2, >0] 340 ELSE 341 zfmdt = rhoic * zh_i(ji,jk) !!! internal melting 342 zdE = 0._wp ! All the layer melts if t_i(jk) = Tmelt (i.e. zdE = 0) 343 END IF 346 344 347 345 zdeltah(ji,jk) = - zfmdt * r1_rhoic ! Melt of layer jk [m, <0] … … 358 356 359 357 ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 360 sfx_sum_1d(ji) = sfx_sum_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic* r1_rdtice358 sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 361 359 362 360 ! Contribution to heat flux [W.m-2], < 0 … … 405 403 ! -> need for an iterative procedure, which converges quickly 406 404 407 IF ( nn_icesal == 2 ) THEN 408 num_iter_max = 5 409 ELSE 410 num_iter_max = 1 411 ENDIF 412 413 ! Just to be sure that enthalpy at nlay_i+1 is null 414 DO ji = kideb, kiut 415 q_i_1d(ji,nlay_i+1) = 0._wp 416 END DO 405 num_iter_max = 1 406 IF( nn_icesal == 2 ) num_iter_max = 5 417 407 418 408 ! Iterative procedure … … 483 473 484 474 ! Contribution to salt flux, <0 485 sfx_bog_1d(ji) = sfx_bog_1d(ji) + s_i_new(ji) * a_i_1d(ji) * zfmdt* r1_rdtice475 sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * s_i_new(ji) * r1_rdtice 486 476 487 477 ! Contribution to mass flux, <0 … … 500 490 DO jk = nlay_i, 1, -1 501 491 DO ji = kideb, kiut 502 IF( zf_tt(ji) > =0._wp .AND. jk > icount(ji) ) THEN ! do not calculate where layer has already disappeared by surface melting492 IF( zf_tt(ji) > 0._wp .AND. jk > icount(ji) ) THEN ! do not calculate where layer has already disappeared by surface melting 503 493 504 494 ztmelts = - tmut * s_i_1d(ji,jk) + rt0 ! Melting point of layer jk (K) … … 524 514 525 515 ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 526 sfx_res_1d(ji) = sfx_res_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic* r1_rdtice516 sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 527 517 528 518 ! Contribution to mass flux … … 559 549 560 550 ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 561 sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic* r1_rdtice551 sfx_bom_1d(ji) = sfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 562 552 563 553 ! Total heat flux used in this process [W.m-2], >0 … … 595 585 zdeltah (ji,1) = - rswitch * zq_rema(ji) / MAX( q_s_1d(ji,1), epsi20 ) 596 586 zdeltah (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_1d(ji) ) ) ! bound melting 597 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,1)598 587 dh_s_tot (ji) = dh_s_tot(ji) + zdeltah(ji,1) 599 588 ht_s_1d (ji) = ht_s_1d(ji) + zdeltah(ji,1) … … 622 611 dh_snowice(ji) = MAX( 0._wp , ( rhosn * ht_s_1d(ji) + (rhoic-rau0) * ht_i_1d(ji) ) / ( rhosn+rau0-rhoic ) ) 623 612 624 ht_i_1d(ji) 625 ht_s_1d(ji) 613 ht_i_1d(ji) = ht_i_1d(ji) + dh_snowice(ji) 614 ht_s_1d(ji) = ht_s_1d(ji) - dh_snowice(ji) 626 615 627 616 ! Salinity of snow ice … … 669 658 ! Update temperature, energy 670 659 !------------------------------------------- 671 !clem bug: we should take snow into account here672 660 DO ji = kideb, kiut 673 661 rswitch = 1.0 - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) ) … … 688 676 WHERE( ht_i_1d == 0._wp ) a_i_1d = 0._wp 689 677 690 CALL wrk_dealloc( jpij, z h_s, zqprec, zq_su, zq_bo, zf_tt, zq_rema )678 CALL wrk_dealloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema ) 691 679 CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s ) 692 680 CALL wrk_dealloc( jpij, nlay_i+1, zdeltah, zh_i ) -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
r5146 r5167 170 170 CALL wrk_alloc( jpij, isnow, ztsub, ztsubit, zh_i, zh_s, zfsw ) 171 171 CALL wrk_alloc( jpij, zf, dzf, zqns_ice_b, zerrit, zdifcase, zftrice, zihic, zghe ) 172 CALL wrk_alloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart=0)173 CALL wrk_alloc( jpij, nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart=0)174 CALL wrk_alloc( jpij, 175 CALL wrk_alloc( jpij, nlay_i+3,3, ztrid )172 CALL wrk_alloc( jpij,nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart=0 ) 173 CALL wrk_alloc( jpij,nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart=0 ) 174 CALL wrk_alloc( jpij,nlay_i+3, zindterm, zindtbis, zdiagbis ) 175 CALL wrk_alloc( jpij,nlay_i+3,3, ztrid ) 176 176 177 177 CALL wrk_alloc( jpij, zdq, zq_ini, zhfx_err ) … … 772 772 ENDIF 773 773 hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err(ji) * a_i_1d(ji) 774 775 ! total heat that is sent to the ocean (i.e. not used in the heat diffusion equation) 776 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err(ji) * a_i_1d(ji) 774 777 END DO 775 776 hfx_err_dif_1d(:) = hfx_err_dif_1d(:) - zhfx_err(:) * a_i_1d(:)777 778 778 779 !----------------------------------------- … … 787 788 & ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji) 788 789 ENDIF 789 END DO 790 791 ! --- compute diagnostic net heat flux at the surface of the snow-ice system (W.m-2) 792 DO ji = kideb, kiut 793 ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 794 hfx_in (ii,ij) = hfx_in (ii,ij) + a_i_1d(ji) * ( qsr_ice_1d(ji) + qns_ice_1d(ji) ) 795 END DO 796 790 ! correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation 791 hfx_dif_1d(ji) = hfx_dif_1d(ji) - zhfx_err(ji) * a_i_1d(ji) 792 END DO 797 793 ! 798 794 CALL wrk_dealloc( jpij, numeqmin, numeqmax ) 799 795 CALL wrk_dealloc( jpij, isnow, ztsub, ztsubit, zh_i, zh_s, zfsw ) 800 796 CALL wrk_dealloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zghe ) 801 CALL wrk_dealloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, & 802 & ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 ) 803 CALL wrk_dealloc( jpij, nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart = 0 ) 804 CALL wrk_dealloc( jpij, nlay_i+3, zindterm, zindtbis, zdiagbis ) 805 CALL wrk_dealloc( jpij, nlay_i+3, 3, ztrid ) 797 CALL wrk_dealloc( jpij,nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 ) 798 CALL wrk_dealloc( jpij,nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart = 0 ) 799 CALL wrk_dealloc( jpij,nlay_i+3, zindterm, zindtbis, zdiagbis ) 800 CALL wrk_dealloc( jpij,nlay_i+3,3, ztrid ) 806 801 CALL wrk_dealloc( jpij, zdq, zq_ini, zhfx_err ) 807 802 -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
r5134 r5167 31 31 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 32 32 USE limthd_ent 33 USE limvar 33 34 34 35 IMPLICIT NONE … … 122 123 CALL wrk_alloc( jpi,jpj, zvrel ) 123 124 125 CALL lim_var_agg(1) 126 CALL lim_var_glo2eqv 124 127 !------------------------------------------------------------------------------| 125 128 ! 2) Convert units for ice internal energy -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
r5138 r5167 112 112 113 113 !--- Thickness correction init. ------------------------------- 114 CALL lim_var_glo2eqv115 114 zatold(:,:) = SUM( a_i(:,:,:), dim=3 ) 115 DO jl = 1, jpl 116 DO jj = 1, jpj 117 DO ji = 1, jpi 118 rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) ) 119 ht_i (ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 120 ht_s (ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 121 END DO 122 END DO 123 END DO 116 124 !--------------------------------------------------------------------- 117 ! Record max of the surrounding ice thicknesses for correction in limupdate125 ! Record max of the surrounding ice thicknesses for correction 118 126 ! in case advection creates ice too thick. 119 127 !--------------------------------------------------------------------- … … 147 155 CALL ctl_warn( 'lim_trp: ncfl= ', TRIM(cltmp), 'advective ice time-step using a split in sub-time-step ') 148 156 ELSE 149 157 ! WRITE(numout,*) 'lim_trp : CFL criterion for ice advection is always smaller than 1/2 ' 150 158 ENDIF 151 159 ENDIF … … 228 236 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl), & 229 237 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 230 231 238 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0oi (:,:,jl), sxage(:,:,jl), & !--- ice age --- 232 239 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) … … 345 352 !!gm & cr 346 353 354 ! --- diags --- 355 DO jj = 1, jpj 356 DO ji = 1, jpi 357 diag_trp_ei(ji,jj) = ( SUM( e_i(ji,jj,1:nlay_i,:) ) - zeiold(ji,jj) ) * r1_rdtice 358 diag_trp_es(ji,jj) = ( SUM( e_s(ji,jj,1:nlay_s,:) ) - zesold(ji,jj) ) * r1_rdtice 359 360 diag_trp_vi (ji,jj) = SUM( v_i(ji,jj,:) - zviold(ji,jj,:) ) * r1_rdtice 361 diag_trp_vs (ji,jj) = SUM( v_s(ji,jj,:) - zvsold(ji,jj,:) ) * r1_rdtice 362 diag_trp_smv(ji,jj) = SUM( smv_i(ji,jj,:) - zsmvold(ji,jj,:) ) * r1_rdtice 363 END DO 364 END DO 365 347 366 ! zap small areas 348 367 CALL lim_var_zapsmall 349 368 350 369 !--- Thickness correction in case too high -------------------------------------------------------- 351 CALL lim_var_glo2eqv352 370 DO jl = 1, jpl 353 371 DO jj = 1, jpj … … 356 374 IF ( v_i(ji,jj,jl) > 0._wp ) THEN 357 375 376 rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) ) 377 ht_i (ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 378 ht_s (ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 379 358 380 zvi = v_i (ji,jj,jl) 359 381 zvs = v_s (ji,jj,jl) … … 365 387 366 388 IF ( ( zdv > 0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) .AND. zatold(ji,jj) < 0.80 ) .OR. & 367 & ( zdv <= 0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) ) ) THEN 389 & ( zdv <= 0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) ) ) THEN 368 390 369 391 rswitch = MAX( 0._wp, SIGN( 1._wp, zhimax(ji,jj,jl) - epsi20 ) ) … … 405 427 ENDIF 406 428 407 ! --- diags ---408 DO jj = 1, jpj409 DO ji = 1, jpi410 diag_trp_ei(ji,jj) = ( SUM( e_i(ji,jj,1:nlay_i,:) ) - zeiold(ji,jj) ) * r1_rdtice411 diag_trp_es(ji,jj) = ( SUM( e_s(ji,jj,1:nlay_s,:) ) - zesold(ji,jj) ) * r1_rdtice412 413 diag_trp_vi (ji,jj) = SUM( v_i(ji,jj,:) - zviold(ji,jj,:) ) * r1_rdtice414 diag_trp_vs (ji,jj) = SUM( v_s(ji,jj,:) - zvsold(ji,jj,:) ) * r1_rdtice415 diag_trp_smv(ji,jj) = SUM( smv_i(ji,jj,:) - zsmvold(ji,jj,:) ) * r1_rdtice416 END DO417 END DO418 419 429 ! --- agglomerate variables ----------------- 420 430 vt_i (:,:) = 0._wp … … 443 453 444 454 ENDIF 445 446 CALL lim_var_glo2eqv ! equivalent variables, requested for rafting447 455 448 456 ! ------------------------------------------------- -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90
r5134 r5167 69 69 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate1', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 70 70 71 CALL lim_var_glo2eqv72 71 !---------------------------------------------------- 73 72 ! ice concentration should not exceed amax … … 88 87 END DO 89 88 90 !----------------------------------------------------91 ! Rebin categories with thickness out of bounds92 !----------------------------------------------------93 IF ( jpl > 1 ) CALL lim_itd_th_reb(1, jpl)94 95 !-----------------96 ! zap small values97 !-----------------98 CALL lim_var_zapsmall99 100 89 !--------------------- 101 90 ! Ice salinity bounds … … 106 95 DO ji = 1, jpi 107 96 zsal = smv_i(ji,jj,jl) 108 smv_i(ji,jj,jl) = sm_i(ji,jj,jl) * v_i(ji,jj,jl)109 97 ! salinity stays in bounds 110 98 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) ) … … 117 105 ENDIF 118 106 107 !---------------------------------------------------- 108 ! Rebin categories with thickness out of bounds 109 !---------------------------------------------------- 110 IF ( jpl > 1 ) CALL lim_itd_th_reb(1, jpl) 111 112 !----------------- 113 ! zap small values 114 !----------------- 115 CALL lim_var_zapsmall 116 117 ! ------------------------------------------------- 118 ! Diagnostics 119 ! ------------------------------------------------- 120 DO jl = 1, jpl 121 afx_dyn(:,:) = afx_dyn(:,:) + ( a_i(:,:,jl) - a_i_b(:,:,jl) ) * r1_rdtice 122 END DO 123 124 DO jj = 1, jpj 125 DO ji = 1, jpi 126 ! heat content variation (W.m-2) 127 diag_heat(ji,jj) = - ( SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) ) + & 128 & SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) ) & 129 & ) * r1_rdtice 130 ! salt, volume 131 diag_smvi(ji,jj) = SUM( smv_i(ji,jj,:) - smv_i_b(ji,jj,:) ) * rhoic * r1_rdtice 132 diag_vice(ji,jj) = SUM( v_i (ji,jj,:) - v_i_b (ji,jj,:) ) * rhoic * r1_rdtice 133 diag_vsnw(ji,jj) = SUM( v_s (ji,jj,:) - v_s_b (ji,jj,:) ) * rhosn * r1_rdtice 134 END DO 135 END DO 136 119 137 ! conservation test 120 138 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limupdate1', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 121 122 ! -------------------------------------------------123 ! Diagnostics124 ! -------------------------------------------------125 DO jl = 1, jpl126 afx_dyn(:,:) = afx_dyn(:,:) + ( a_i(:,:,jl) - a_i_b(:,:,jl) ) * r1_rdtice127 END DO128 129 ! heat content variation (W.m-2)130 DO jj = 1, jpj131 DO ji = 1, jpi132 diag_heat_dhc(ji,jj) = - ( SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) ) + &133 & SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) ) &134 & ) * r1_rdtice135 END DO136 END DO137 139 138 140 ! ------------------------------------------------- -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90
r5134 r5167 72 72 ! Constrain the thickness of the smallest category above himin 73 73 !---------------------------------------------------------------------- 74 CALL lim_var_glo2eqv75 74 DO jj = 1, jpj 76 75 DO ji = 1, jpi 76 rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,1) - epsi20 ) ) !0 if no ice and 1 if yes 77 ht_i(ji,jj,1) = v_i (ji,jj,1) / MAX( a_i(ji,jj,1) , epsi20 ) * rswitch 77 78 IF( v_i(ji,jj,1) > 0._wp .AND. ht_i(ji,jj,1) < rn_himin ) THEN 78 79 a_i (ji,jj,1) = a_i(ji,jj,1) * ht_i(ji,jj,1) / rn_himin … … 98 99 END DO 99 100 END DO 100 101 !----------------------------------------------------102 ! Rebin categories with thickness out of bounds103 !----------------------------------------------------104 IF ( jpl > 1 ) CALL lim_itd_th_reb( 1, jpl )105 106 !-----------------107 ! zap small values108 !-----------------109 CALL lim_var_zapsmall110 101 111 102 !--------------------- … … 117 108 DO ji = 1, jpi 118 109 zsal = smv_i(ji,jj,jl) 119 smv_i(ji,jj,jl) = sm_i(ji,jj,jl) * v_i(ji,jj,jl)120 110 ! salinity stays in bounds 121 111 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) ) … … 127 117 END DO 128 118 ENDIF 119 120 !---------------------------------------------------- 121 ! Rebin categories with thickness out of bounds 122 !---------------------------------------------------- 123 IF ( jpl > 1 ) CALL lim_itd_th_reb( 1, jpl ) 124 125 !----------------- 126 ! zap small values 127 !----------------- 128 CALL lim_var_zapsmall 129 129 130 130 !------------------------------------------------------------------------------ … … 150 150 v_ice(:,:) = v_ice(:,:) * vmask(:,:,1) 151 151 152 ! for outputs153 CALL lim_var_glo2eqv ! equivalent variables (outputs)154 CALL lim_var_agg(2) ! aggregate ice thickness categories155 156 ! conservation test157 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate2', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)158 159 152 ! ------------------------------------------------- 160 153 ! Diagnostics … … 165 158 afx_tot = afx_thd + afx_dyn 166 159 167 ! heat content variation (W.m-2)168 160 DO jj = 1, jpj 169 161 DO ji = 1, jpi 170 diag_heat_dhc(ji,jj) = diag_heat_dhc(ji,jj) - & 171 & ( SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) ) + & 172 & SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) ) & 173 & ) * r1_rdtice 174 END DO 175 END DO 162 ! heat content variation (W.m-2) 163 diag_heat(ji,jj) = diag_heat(ji,jj) - & 164 & ( SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) ) + & 165 & SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) ) & 166 & ) * r1_rdtice 167 ! salt, volume 168 diag_smvi(ji,jj) = diag_smvi(ji,jj) + SUM( smv_i(ji,jj,:) - smv_i_b(ji,jj,:) ) * rhoic * r1_rdtice 169 diag_vice(ji,jj) = diag_vice(ji,jj) + SUM( v_i (ji,jj,:) - v_i_b (ji,jj,:) ) * rhoic * r1_rdtice 170 diag_vsnw(ji,jj) = diag_vsnw(ji,jj) + SUM( v_s (ji,jj,:) - v_s_b (ji,jj,:) ) * rhosn * r1_rdtice 171 END DO 172 END DO 173 174 ! conservation test 175 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limupdate2', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 176 177 ! for outputs 178 CALL lim_var_glo2eqv 179 CALL lim_var_agg(2) 176 180 177 181 ! ------------------------------------------------- -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
r5134 r5167 124 124 DO ji = 1, jpi 125 125 et_s(ji,jj) = et_s(ji,jj) + e_s(ji,jj,1,jl) ! snow heat content 126 rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi 10 ) )127 smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi 10 ) * rswitch ! ice salinity128 rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi 10 ) )129 ot_i(ji,jj) = ot_i(ji,jj) + oa_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi 10 ) * rswitch ! ice age126 rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi20 ) ) 127 smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi20 ) * rswitch ! ice salinity 128 rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi20 ) ) 129 ot_i(ji,jj) = ot_i(ji,jj) + oa_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi20 ) * rswitch ! ice age 130 130 END DO 131 131 END DO … … 161 161 DO jj = 1, jpj 162 162 DO ji = 1, jpi 163 rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi 10 ) ) !0 if no ice and 1 if yes164 ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi 10 ) * rswitch165 ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi 10 ) * rswitch166 o_i(ji,jj,jl) = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi 10 ) * rswitch163 rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) ) !0 if no ice and 1 if yes 164 ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 165 ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 166 o_i(ji,jj,jl) = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 167 167 END DO 168 168 END DO … … 173 173 DO jj = 1, jpj 174 174 DO ji = 1, jpi 175 rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi 10 ) ) !0 if no ice and 1 if yes176 sm_i(ji,jj,jl) = smv_i(ji,jj,jl) / MAX( v_i(ji,jj,jl) , epsi 10 ) * rswitch175 rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi20 ) ) !0 if no ice and 1 if yes 176 sm_i(ji,jj,jl) = smv_i(ji,jj,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * rswitch 177 177 END DO 178 178 END DO … … 228 228 ! Mean temperature 229 229 !------------------- 230 vt_i (:,:) = 0._wp 231 DO jl = 1, jpl 232 vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl) 233 END DO 234 230 235 tm_i(:,:) = 0._wp 231 236 DO jl = 1, jpl … … 233 238 DO jj = 1, jpj 234 239 DO ji = 1, jpi 235 rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) ) 236 tm_i(ji,jj) = tm_i(ji,jj) + rswitch * t_i(ji,jj,jk,jl) * v_i(ji,jj,jl) & 237 & / ( REAL(nlay_i,wp) * MAX( vt_i(ji,jj) , epsi10 ) ) 240 rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi20 ) ) 241 tm_i(ji,jj) = tm_i(ji,jj) + ( 1._wp - rswitch ) * rt0 + & 242 & rswitch * t_i(ji,jj,jk,jl) * v_i(ji,jj,jl) & 243 & / ( REAL(nlay_i,wp) * MAX( vt_i(ji,jj) , epsi20 ) ) 238 244 END DO 239 245 END DO … … 305 311 DO jj = 1, jpj 306 312 DO ji = 1, jpi 307 z_slope_s(ji,jj,jl) = 2._wp * sm_i(ji,jj,jl) / MAX( epsi10 , ht_i(ji,jj,jl) ) 313 rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i(ji,jj,jl) - epsi20 ) ) 314 z_slope_s(ji,jj,jl) = rswitch * 2._wp * sm_i(ji,jj,jl) / MAX( epsi20 , ht_i(ji,jj,jl) ) 308 315 END DO 309 316 END DO … … 379 386 380 387 ! Mean sea ice temperature 388 vt_i (:,:) = 0._wp 389 DO jl = 1, jpl 390 vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl) 391 END DO 392 381 393 tm_i(:,:) = 0._wp 382 394 DO jl = 1, jpl … … 384 396 DO jj = 1, jpj 385 397 DO ji = 1, jpi 386 rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi 10) )398 rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi06 ) ) 387 399 tm_i(ji,jj) = tm_i(ji,jj) + rswitch * t_i(ji,jj,jk,jl) * v_i(ji,jj,jl) & 388 & * r1_nlay_i / MAX( vt_i(ji,jj) , epsi 10)400 & * r1_nlay_i / MAX( vt_i(ji,jj) , epsi06 ) 389 401 END DO 390 402 END DO … … 409 421 !!------------------------------------------------------------------ 410 422 ! 423 vt_i (:,:) = 0._wp 424 DO jl = 1, jpl 425 vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl) 426 END DO 427 411 428 bv_i(:,:) = 0._wp 412 429 DO jl = 1, jpl … … 417 434 zbvi = - rswitch * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rt0, - epsi10 ) & 418 435 & * v_i(ji,jj,jl) * r1_nlay_i 419 rswitch = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi 10 ) ) )420 bv_i(ji,jj) = bv_i(ji,jj) + rswitch * zbvi / MAX( vt_i(ji,jj) , epsi 10 )436 rswitch = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi20 ) ) ) 437 bv_i(ji,jj) = bv_i(ji,jj) + rswitch * zbvi / MAX( vt_i(ji,jj) , epsi20 ) 421 438 END DO 422 439 END DO … … 460 477 ! 461 478 DO ji = kideb, kiut ! Slope of the linear profile zs_zero 462 z_slope_s(ji) = 2._wp * sm_i_1d(ji) / MAX( epsi10 , ht_i_1d(ji) ) 479 rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi20 ) ) 480 z_slope_s(ji) = rswitch * 2._wp * sm_i_1d(ji) / MAX( epsi20 , ht_i_1d(ji) ) 463 481 END DO 464 482 -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90
r5123 r5167 232 232 CALL iom_put ('hfxopw' , hfx_opw(:,:) ) ! 233 233 CALL iom_put ('hfxtur' , fhtur(:,:) * at_i(:,:) ) ! turbulent heat flux at ice base 234 CALL iom_put ('hfxdhc' , diag_heat _dhc(:,:)) ! Heat content variation in snow and ice234 CALL iom_put ('hfxdhc' , diag_heat(:,:) ) ! Heat content variation in snow and ice 235 235 CALL iom_put ('hfxspr' , hfx_spr(:,:) ) ! Heat content of snow precip 236 236 -
trunk/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90
r5146 r5167 20 20 ! !!! ** ice-thermo namelist (namicethd) ** 21 21 REAL(wp), PUBLIC :: rn_himin !: minimum ice thickness 22 REAL(wp), PUBLIC :: parsub !: switch for snow sublimation or not23 22 REAL(wp), PUBLIC :: rn_maxfrazb !: maximum portion of frazil ice collecting at the ice bottom 24 23 REAL(wp), PUBLIC :: rn_vfrazb !: threshold drift speed for collection of bottom frazil ice … … 140 139 !!---------------------------------------------------------------------! 141 140 142 ALLOCATE( npb (jpij) , nplm (jpij) , npac (jpij), & 143 ! ! 144 & qlead_1d (jpij) , ftr_ice_1d (jpij) , & 145 & qsr_ice_1d (jpij) , & 146 & fr1_i0_1d(jpij) , fr2_i0_1d(jpij) , qns_ice_1d(jpij) , & 141 ALLOCATE( npb (jpij) , nplm (jpij) , npac (jpij) , & 142 & qlead_1d (jpij) , ftr_ice_1d(jpij) , qsr_ice_1d (jpij) , & 143 & fr1_i0_1d(jpij) , fr2_i0_1d (jpij) , qns_ice_1d(jpij) , & 147 144 & t_bo_1d (jpij) , & 148 145 & hfx_sum_1d(jpij) , hfx_bom_1d(jpij) ,hfx_bog_1d(jpij) , & … … 152 149 & hfx_res_1d(jpij) , hfx_err_rem_1d(jpij) , hfx_err_dif_1d(jpij) , STAT=ierr(1) ) 153 150 ! 154 ALLOCATE( sprecip_1d (jpij) , frld_1d (jpij) , at_i_1d (jpij) , & 155 & fhtur_1d (jpij) , wfx_snw_1d (jpij) , wfx_spr_1d (jpij) , & 156 & fhld_1d (jpij) , wfx_sub_1d (jpij) , wfx_bog_1d(jpij) , wfx_bom_1d(jpij) , & 157 & wfx_sum_1d(jpij) , wfx_sni_1d (jpij) , wfx_opw_1d (jpij) , wfx_res_1d (jpij) , & 158 & dqns_ice_1d(jpij) , qla_ice_1d (jpij) , dqla_ice_1d(jpij) , & 159 & tatm_ice_1d(jpij) , & 160 & i0 (jpij) , & 161 & sfx_bri_1d (jpij) , sfx_bog_1d (jpij) , sfx_bom_1d (jpij) ,sfx_sum_1d (jpij) , & 162 & sfx_sni_1d (jpij) , sfx_opw_1d (jpij) , sfx_res_1d (jpij) , & 163 & dsm_i_fl_1d(jpij) , dsm_i_gd_1d(jpij) , dsm_i_se_1d(jpij) , & 151 ALLOCATE( sprecip_1d (jpij) , frld_1d (jpij) , at_i_1d (jpij) , & 152 & fhtur_1d (jpij) , wfx_snw_1d (jpij) , wfx_spr_1d (jpij) , & 153 & fhld_1d (jpij) , wfx_sub_1d (jpij) , wfx_bog_1d (jpij) , wfx_bom_1d(jpij) , & 154 & wfx_sum_1d(jpij) , wfx_sni_1d (jpij) , wfx_opw_1d (jpij) , wfx_res_1d(jpij) , & 155 & dqns_ice_1d(jpij) , qla_ice_1d (jpij) , dqla_ice_1d(jpij) , & 156 & tatm_ice_1d(jpij) , i0 (jpij) , & 157 & sfx_bri_1d (jpij) , sfx_bog_1d (jpij) , sfx_bom_1d (jpij) , sfx_sum_1d (jpij), & 158 & sfx_sni_1d (jpij) , sfx_opw_1d (jpij) , sfx_res_1d (jpij) , & 159 & dsm_i_fl_1d(jpij) , dsm_i_gd_1d(jpij) , dsm_i_se_1d(jpij) , & 164 160 & dsm_i_si_1d(jpij) , hicol_1d (jpij) , STAT=ierr(2) ) 165 161 ! 166 ALLOCATE( t_su_1d (jpij) , a_i_1d (jpij) , ht_i_1d(jpij) , &167 & ht_s_1d 162 ALLOCATE( t_su_1d (jpij) , a_i_1d (jpij) , ht_i_1d (jpij) , & 163 & ht_s_1d (jpij) , fc_su (jpij) , fc_bo_i (jpij) , & 168 164 & dh_s_tot (jpij) , dh_i_surf(jpij) , dh_i_bott(jpij) , & 169 & dh_snowice(jpij) , & 170 & sm_i_1d (jpij) , s_i_new (jpij) , & 171 & t_s_1d(jpij,nlay_s), & 172 & t_i_1d(jpij,nlay_i+1), s_i_1d(jpij,nlay_i+1) , & 173 & q_i_1d(jpij,nlay_i+1), q_s_1d(jpij,nlay_i+1) , & 165 & dh_snowice(jpij) , sm_i_1d (jpij) , s_i_new (jpij) , & 166 & t_s_1d(jpij,nlay_s) , t_i_1d(jpij,nlay_i) , s_i_1d(jpij,nlay_i) , & 167 & q_i_1d(jpij,nlay_i+1) , q_s_1d(jpij,nlay_s) , & 174 168 & qh_i_old(jpij,0:nlay_i+1), h_i_old(jpij,0:nlay_i+1) , STAT=ierr(3)) 175 169 ! -
trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim.F90
r5143 r5167 28 28 USE ice ! LIM_3 ice variables 29 29 USE dom_ice ! sea-ice domain 30 USE limvar 30 31 #endif 31 32 USE par_oce ! ocean parameters … … 61 62 INTEGER :: ib_bdy ! Loop index 62 63 64 #if defined key_lim3 65 CALL lim_var_glo2eqv 66 #endif 67 63 68 DO ib_bdy=1, nb_bdy 64 69 … … 73 78 74 79 END DO 80 81 #if defined key_lim3 82 CALL lim_var_zapsmall 83 CALL lim_var_agg(1) 84 #endif 75 85 76 86 END SUBROUTINE bdy_ice_lim -
trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
r5146 r5167 184 184 numit = numit + nn_fsbc ! Ice model time step 185 185 ! 186 CALL sbc_lim_ update! Store previous ice values186 CALL sbc_lim_bef ! Store previous ice values 187 187 188 188 CALL sbc_lim_diag0 ! set diag of mass, heat and salt fluxes to 0 … … 202 202 203 203 #if defined key_bdy 204 CALL lim_var_glo2eqv205 204 CALL bdy_ice_lim( kt ) ! bdy ice thermo 206 CALL lim_var_zapsmall207 CALL lim_var_agg(1)208 205 IF( ln_icectl ) CALL lim_prt( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' ) 209 206 #endif … … 212 209 ENDIF 213 210 214 CALL sbc_lim_ update! Store previous ice values211 CALL sbc_lim_bef ! Store previous ice values 215 212 216 213 ! ---------------------------------------------- 217 214 ! ice thermodynamics 218 215 ! ---------------------------------------------- 219 CALL lim_var_glo2eqv220 216 CALL lim_var_agg(1) 221 217 … … 248 244 ! 249 245 IF( lrst_ice ) CALL lim_rst_write( kt ) ! Ice restart file 250 CALL lim_var_glo2eqv ! ??? 251 ! 252 IF( ln_icectl ) CALL lim_ctl( kt ) ! alerts in case of model crash 246 ! 247 IF( ln_icectl ) CALL lim_ctl( kt ) ! alerts in case of model crash 253 248 ! 254 249 CALL wrk_dealloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice ) … … 389 384 r1_nlay_s = 1._wp / REAL( nlay_s, wp ) 390 385 ! 386 #if defined key_bdy 387 IF( lwp .AND. ln_limdiahsb ) CALL ctl_warn('online conservation check activated but it does not work with BDY') 388 #endif 389 ! 391 390 END SUBROUTINE ice_run 392 391 … … 555 554 END SUBROUTINE ice_lim_flx 556 555 557 SUBROUTINE sbc_lim_ update558 !!---------------------------------------------------------------------- 559 !! *** ROUTINE sbc_lim_ update***556 SUBROUTINE sbc_lim_bef 557 !!---------------------------------------------------------------------- 558 !! *** ROUTINE sbc_lim_bef *** 560 559 !! 561 560 !! ** purpose : store ice variables at "before" time step … … 571 570 v_ice_b(:,:) = v_ice(:,:) 572 571 573 END SUBROUTINE sbc_lim_ update572 END SUBROUTINE sbc_lim_bef 574 573 575 574 SUBROUTINE sbc_lim_diag0 … … 602 601 hfx_spr(:,:) = 0._wp ; hfx_dif(:,:) = 0._wp 603 602 hfx_err(:,:) = 0._wp ; hfx_err_rem(:,:) = 0._wp 604 hfx_err_dif(:,:) = 0._wp 603 hfx_err_dif(:,:) = 0._wp ; 605 604 606 605 afx_tot(:,:) = 0._wp ; 607 606 afx_dyn(:,:) = 0._wp ; afx_thd(:,:) = 0._wp 608 607 609 diag_heat_dhc(:,:) = 0._wp ; 608 diag_heat(:,:) = 0._wp ; diag_smvi(:,:) = 0._wp ; 609 diag_vice(:,:) = 0._wp ; diag_vsnw(:,:) = 0._wp ; 610 610 611 611 END SUBROUTINE sbc_lim_diag0 … … 636 636 637 637 fice_ice_ave (:,:) = 0.0_wp 638 WHERE ( at_i (:,:) .GT.0.0_wp ) fice_ice_ave (:,:) = fice_cell_ave ( ptab (:,:,:)) / at_i (:,:)638 WHERE ( at_i (:,:) > 0.0_wp ) fice_ice_ave (:,:) = fice_cell_ave ( ptab (:,:,:)) / at_i (:,:) 639 639 640 640 END FUNCTION fice_ice_ave
Note: See TracChangeset
for help on using the changeset viewer.
