- Timestamp:
- 2014-11-28T18:24:01+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
r4624 r4924 5 5 !!====================================================================== 6 6 !! History : LIM ! 2006-02 (M. Vancoppenolle) Original code 7 !! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in smsw & sfx_ mec7 !! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in smsw & sfx_dyn 8 8 !! 4.0 ! 2011-02 (G. Madec) dynamical allocation 9 9 !!---------------------------------------------------------------------- … … 22 22 USE limthd_lac ! LIM 23 23 USE limvar ! LIM 24 USE limcons ! LIM25 24 USE in_out_manager ! I/O manager 26 25 USE lbclnk ! lateral boundary condition - MPP exchanges … … 30 29 ! Check budget (Rousset) 31 30 USE iom ! I/O manager 32 USE lib_fortran ! glob_sum31 USE lib_fortran ! glob_sum 33 32 USE limdiahsb 34 USE timing ! Timing 33 USE timing ! Timing 34 USE limcons ! conservation tests 35 35 36 36 IMPLICIT NONE … … 143 143 REAL(wp), POINTER, DIMENSION(:,:) :: esnow_mlt ! energy needed to melt snow in ocean (J m-2) 144 144 REAL(wp), POINTER, DIMENSION(:,:) :: vt_i_init, vt_i_final ! ice volume summed over categories 145 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) 146 REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 147 ! mass and salt flux (clem) 148 REAL(wp), POINTER, DIMENSION(:,:,:) :: zviold, zvsold, zsmvold ! old ice volume... 145 ! 146 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 149 147 !!----------------------------------------------------------------------------- 150 148 IF( nn_timing == 1 ) CALL timing_start('limitd_me') 151 149 152 150 CALL wrk_alloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final ) 153 154 CALL wrk_alloc( jpi, jpj, jpl, zviold, zvsold, zsmvold ) ! clem155 156 IF( numit == nstart ) CALL lim_itd_me_init ! Initialization (first time-step only)157 151 158 152 IF(ln_ctl) THEN … … 162 156 163 157 IF( ln_limdyn ) THEN ! Start ridging and rafting ! 164 ! ------------------------------- 165 !- check conservation (C Rousset) 166 IF (ln_limdiahsb) THEN 167 zchk_v_i_b = glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 168 zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 169 zchk_fw_b = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) 170 zchk_fs_b = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) 171 ENDIF 172 !- check conservation (C Rousset) 173 ! ------------------------------- 174 175 ! mass and salt flux init (clem) 176 zviold(:,:,:) = v_i(:,:,:) 177 zvsold(:,:,:) = v_s(:,:,:) 178 zsmvold(:,:,:) = smv_i(:,:,:) 158 159 ! conservation test 160 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 179 161 180 162 !-----------------------------------------------------------------------------! … … 362 344 ! 5) Heat, salt and freshwater fluxes 363 345 !-----------------------------------------------------------------------------! 364 fmmec(ji,jj) = fmmec(ji,jj) + msnow_mlt(ji,jj) * r1_rdtice ! fresh water source for ocean365 fhmec(ji,jj) = fhmec(ji,jj) + esnow_mlt(ji,jj) * r1_rdtice ! heat sink for ocean346 wfx_snw(ji,jj) = wfx_snw(ji,jj) + msnow_mlt(ji,jj) * r1_rdtice ! fresh water source for ocean 347 hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + esnow_mlt(ji,jj) * unit_fac / area(ji,jj) * r1_rdtice ! heat sink for ocean (<0, W.m-2) 366 348 367 349 END DO … … 399 381 CALL lim_itd_me_zapsmall 400 382 401 !--------------------------------402 ! Update mass/salt fluxes (clem)403 !--------------------------------404 DO jl = 1, jpl405 DO jj = 1, jpj406 DO ji = 1, jpi407 diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * r1_rdtice408 rdm_ice(ji,jj) = rdm_ice(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * rhoic409 rdm_snw(ji,jj) = rdm_snw(ji,jj) + ( v_s(ji,jj,jl) - zvsold(ji,jj,jl) ) * rhosn410 sfx_mec(ji,jj) = sfx_mec(ji,jj) - ( smv_i(ji,jj,jl) - zsmvold(ji,jj,jl) ) * rhoic * r1_rdtice411 END DO412 END DO413 END DO414 383 415 384 IF(ln_ctl) THEN ! Control print … … 445 414 ENDIF 446 415 447 ! ------------------------------- 448 !- check conservation (C Rousset) 449 IF (ln_limdiahsb) THEN 450 zchk_fs = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 451 zchk_fw = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b 452 453 zchk_v_i = ( glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) * r1_rdtice 454 zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic ) 455 456 zchk_vmin = glob_min(v_i) 457 zchk_amax = glob_max(SUM(a_i,dim=3)) 458 zchk_amin = glob_min(a_i) 459 460 IF(lwp) THEN 461 IF ( ABS( zchk_v_i ) > 1.e-5 ) WRITE(numout,*) 'violation volume [m3/day] (limitd_me) = ',(zchk_v_i * rday) 462 IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limitd_me) = ',(zchk_smv * rday) 463 IF ( zchk_vmin < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limitd_me) = ',(zchk_vmin * 1.e-3) 464 IF ( zchk_amax > kamax+epsi10 ) WRITE(numout,*) 'violation a_i>amax (limitd_me) = ',zchk_amax 465 IF ( zchk_amin < 0. ) WRITE(numout,*) 'violation a_i<0 (limitd_me) = ',zchk_amin 466 ENDIF 467 ENDIF 468 !- check conservation (C Rousset) 469 ! ------------------------------- 416 ! conservation test 417 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 470 418 471 419 ENDIF ! ln_limdyn=.true. 472 420 ! 473 421 CALL wrk_dealloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final ) 474 !475 CALL wrk_dealloc( jpi, jpj, jpl, zviold, zvsold, zsmvold ) ! clem476 422 ! 477 423 IF( nn_timing == 1 ) CALL timing_stop('limitd_me') … … 670 616 !!---------------------------------------------------------------------! 671 617 INTEGER :: ji,jj, jl ! dummy loop indices 672 INTEGER :: krdg_index !673 618 REAL(wp) :: Gstari, astari, hi, hrmean, zdummy ! local scalar 674 619 REAL(wp), POINTER, DIMENSION(:,:) :: zworka ! temporary array used here … … 746 691 !----------------------------------------------------------------- 747 692 748 krdg_index = 1 749 750 IF( krdg_index == 0 ) THEN !--- Linear formulation (Thorndike et al., 1975) 751 DO jl = 0, ice_cat_bounds(1,2) ! only undeformed ice participates 693 IF( partfun_swi == 0 ) THEN !--- Linear formulation (Thorndike et al., 1975) 694 DO jl = 0, jpl 752 695 DO jj = 1, jpj 753 696 DO ji = 1, jpi … … 772 715 Gsum(:,:,jl) = EXP( -Gsum(:,:,jl) * astari ) * zdummy 773 716 END DO !jl 774 DO jl = 0, ice_cat_bounds(1,2)717 DO jl = 0, jpl 775 718 athorn(:,:,jl) = Gsum(:,:,jl-1) - Gsum(:,:,jl) 776 719 END DO 777 720 ! 778 ENDIF ! krdg_index779 780 IF( raft swi == 1 ) THEN ! Ridging and rafting ice participation functions721 ENDIF ! partfun_swi 722 723 IF( raft_swi == 1 ) THEN ! Ridging and rafting ice participation functions 781 724 ! 782 725 DO jl = 1, jpl … … 794 737 END DO ! jl 795 738 796 ELSE ! raft swi = 0739 ELSE ! raft_swi = 0 797 740 ! 798 741 DO jl = 1, jpl … … 802 745 ENDIF 803 746 804 IF ( raft swi == 1 ) THEN747 IF ( raft_swi == 1 ) THEN 805 748 806 749 IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) .GT. epsi10 ) THEN … … 908 851 INTEGER :: ij ! horizontal index, combines i and j loops 909 852 INTEGER :: icells ! number of cells with aicen > puny 910 REAL(wp) :: zindb , zsrdg2! local scalar853 REAL(wp) :: zindb ! local scalar 911 854 REAL(wp) :: hL, hR, farea, zdummy, zdummy0, ztmelts ! left and right limits of integration 855 REAL(wp) :: zsstK ! SST in Kelvin 912 856 913 857 INTEGER , POINTER, DIMENSION(:) :: indxi, indxj ! compressed indices … … 917 861 918 862 REAL(wp), POINTER, DIMENSION(:,:,:) :: aicen_init, vicen_init ! ice area & volume before ridging 919 REAL(wp), POINTER, DIMENSION(:,:,:) :: vsn on_init, esnon_init ! snow volume & energy before ridging863 REAL(wp), POINTER, DIMENSION(:,:,:) :: vsnwn_init, esnwn_init ! snow volume & energy before ridging 920 864 REAL(wp), POINTER, DIMENSION(:,:,:) :: smv_i_init, oa_i_init ! ice salinity & age before ridging 921 865 … … 952 896 CALL wrk_alloc( jpi, jpj, vrdg1, vrdg2, vsw , srdg1, srdg2, smsw ) 953 897 CALL wrk_alloc( jpi, jpj, afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 ) 954 CALL wrk_alloc( jpi, jpj, jpl, aicen_init, vicen_init, vsn on_init, esnon_init, smv_i_init, oa_i_init )955 CALL wrk_alloc( jpi, jpj, jkmax, eirft, erdg1, erdg2, ersw )956 CALL wrk_alloc( jpi, jpj, jkmax, jpl, eicen_init )898 CALL wrk_alloc( jpi, jpj, jpl, aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init ) 899 CALL wrk_alloc( jpi, jpj, nlay_i+1, eirft, erdg1, erdg2, ersw ) 900 CALL wrk_alloc( jpi, jpj, nlay_i+1, jpl, eicen_init ) 957 901 958 902 ! Conservation check … … 1008 952 aicen_init(:,:,jl) = a_i(:,:,jl) 1009 953 vicen_init(:,:,jl) = v_i(:,:,jl) 1010 vsn on_init(:,:,jl) = v_s(:,:,jl)954 vsnwn_init(:,:,jl) = v_s(:,:,jl) 1011 955 ! 1012 956 smv_i_init(:,:,jl) = smv_i(:,:,jl) … … 1014 958 END DO !jl 1015 959 1016 esn on_init(:,:,:) = e_s(:,:,1,:)960 esnwn_init(:,:,:) = e_s(:,:,1,:) 1017 961 1018 962 DO jl = 1, jpl … … 1091 1035 ! / rafting category n1. 1092 1036 !-------------------------------------------------------------------------- 1093 vrdg1(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj) / ( 1._wp + ridge_por )1037 vrdg1(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj) 1094 1038 vrdg2(ji,jj) = vrdg1(ji,jj) * ( 1. + ridge_por ) 1095 1039 vsw (ji,jj) = vrdg1(ji,jj) * ridge_por 1096 1040 1097 vsrdg(ji,jj) = vsn on_init(ji,jj,jl1) * afrac(ji,jj)1098 esrdg(ji,jj) = esn on_init(ji,jj,jl1) * afrac(ji,jj)1099 srdg1(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) / ( 1._wp + ridge_por )1100 srdg2(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) 1041 vsrdg(ji,jj) = vsnwn_init(ji,jj,jl1) * afrac(ji,jj) 1042 esrdg(ji,jj) = esnwn_init(ji,jj,jl1) * afrac(ji,jj) 1043 srdg1(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) 1044 srdg2(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) !! MV HC 2014 this line seems useless 1101 1045 1102 1046 ! rafting volumes, heat contents ... 1103 1047 virft(ji,jj) = vicen_init(ji,jj,jl1) * afrft(ji,jj) 1104 vsrft(ji,jj) = vsn on_init(ji,jj,jl1) * afrft(ji,jj)1105 esrft(ji,jj) = esn on_init(ji,jj,jl1) * afrft(ji,jj)1048 vsrft(ji,jj) = vsnwn_init(ji,jj,jl1) * afrft(ji,jj) 1049 esrft(ji,jj) = esnwn_init(ji,jj,jl1) * afrft(ji,jj) 1106 1050 smrft(ji,jj) = smv_i_init(ji,jj,jl1) * afrft(ji,jj) 1107 1051 … … 1120 1064 ! Salinity 1121 1065 !------------- 1122 smsw(ji,jj) = sss_m(ji,jj) * vsw(ji,jj) * rhoic / rau0 ! salt content of seawater frozen in voids 1123 1124 zsrdg2 = srdg1(ji,jj) + smsw(ji,jj) ! salt content of new ridge 1125 1126 srdg2(ji,jj) = MIN( s_i_max * vrdg2(ji,jj) , zsrdg2 ) ! impose a maximum salinity 1066 smsw(ji,jj) = vsw(ji,jj) * sss_m(ji,jj) ! salt content of seawater frozen in voids !! MV HC2014 1067 srdg2(ji,jj) = srdg1(ji,jj) + smsw(ji,jj) ! salt content of new ridge 1068 1069 !srdg2(ji,jj) = MIN( s_i_max * vrdg2(ji,jj) , zsrdg2 ) ! impose a maximum salinity 1127 1070 1128 ! ! excess of salt is flushed into the ocean 1129 !sfx_mec(ji,jj) = sfx_mec(ji,jj) + ( zsrdg2 - srdg2(ji,jj) ) * rhoic * r1_rdtice 1130 1131 !rdm_ice(ji,jj) = rdm_ice(ji,jj) + vsw(ji,jj) * rhoic ! gurvan: increase in ice volume du to seawater frozen in voids 1071 sfx_dyn(ji,jj) = sfx_dyn(ji,jj) - smsw(ji,jj) * rhoic * r1_rdtice 1072 wfx_dyn(ji,jj) = wfx_dyn(ji,jj) - vsw (ji,jj) * rhoic * r1_rdtice ! gurvan: increase in ice volume du to seawater frozen in voids 1132 1073 1133 1074 !------------------------------------ … … 1158 1099 & + rhosn*vsrft(ji,jj)*(1.0-fsnowrft) 1159 1100 1160 esnow_mlt(ji,jj) = esnow_mlt(ji,jj) + esrdg(ji,jj)*(1.0-fsnowrdg) & !rafting included 1161 & + esrft(ji,jj)*(1.0-fsnowrft) 1101 ! in 1e-9 Joules (same as e_s) 1102 esnow_mlt(ji,jj) = esnow_mlt(ji,jj) - esrdg(ji,jj)*(1.0-fsnowrdg) & !rafting included 1103 & - esrft(ji,jj)*(1.0-fsnowrft) 1162 1104 1163 1105 !----------------------------------------------------------------- … … 1184 1126 jj = indxj(ij) 1185 1127 ! heat content of ridged ice 1186 erdg1(ji,jj,jk) = eicen_init(ji,jj,jk,jl1) * afrac(ji,jj) / ( 1._wp + ridge_por )1128 erdg1(ji,jj,jk) = eicen_init(ji,jj,jk,jl1) * afrac(ji,jj) 1187 1129 eirft(ji,jj,jk) = eicen_init(ji,jj,jk,jl1) * afrft(ji,jj) 1188 1130 e_i (ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) - erdg1(ji,jj,jk) - eirft(ji,jj,jk) 1189 ! sea water heat content 1190 ztmelts = - tmut * sss_m(ji,jj) + rtt 1191 ! heat content per unit volume 1192 zdummy0 = - rcp * ( sst_m(ji,jj) + rt0 - rtt ) * vsw(ji,jj) 1193 1194 ! corrected sea water salinity 1195 zindb = MAX( 0._wp , SIGN( 1._wp , vsw(ji,jj) - epsi20 ) ) 1196 zdummy = zindb * ( srdg1(ji,jj) - srdg2(ji,jj) ) / MAX( ridge_por * vsw(ji,jj), epsi20 ) 1197 1198 ztmelts = - tmut * zdummy + rtt 1199 ersw(ji,jj,jk) = - rcp * ( ztmelts - rtt ) * vsw(ji,jj) 1200 1201 ! heat flux 1202 fheat_mec(ji,jj) = fheat_mec(ji,jj) + ( ersw(ji,jj,jk) - zdummy0 ) * r1_rdtice 1131 1132 1133 ! enthalpy of the trapped seawater (J/m2, >0) 1134 ! clem: if sst>0, then ersw <0 (is that possible?) 1135 zsstK = sst_m(ji,jj) + rt0 1136 ersw(ji,jj,jk) = - rhoic * vsw(ji,jj) * rcp * ( zsstK - rt0 ) / REAL( nlay_i ) 1137 1138 ! heat flux to the ocean 1139 hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ersw(ji,jj,jk) * r1_rdtice ! > 0 [W.m-2] ocean->ice flux 1203 1140 1204 1141 ! Correct dimensions to avoid big values 1205 ersw(ji,jj,jk) = ersw(ji,jj,jk) * 1.e-09 1206 1207 ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J 1208 ersw (ji,jj,jk) = ersw(ji,jj,jk) * area(ji,jj) * vsw(ji,jj) / REAL( nlay_i ) 1142 ersw(ji,jj,jk) = ersw(ji,jj,jk) / unit_fac 1143 1144 ! Mutliply by ice volume, and divide by number of layers to get heat content in 1.e9 J 1145 ! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean 1146 !! MV HC 2014 1147 ersw (ji,jj,jk) = ersw(ji,jj,jk) * area(ji,jj) 1209 1148 1210 1149 erdg2(ji,jj,jk) = erdg1(ji,jj,jk) + ersw(ji,jj,jk) 1150 1211 1151 END DO ! ij 1212 1152 END DO !jk … … 1253 1193 !------------------------------------------------------------------------------- 1254 1194 ! jl1 looping 1-jpl 1255 DO jl2 = ice_cat_bounds(1,1), ice_cat_bounds(1,2)1195 DO jl2 = 1, jpl 1256 1196 ! over categories to which ridged ice is transferred 1257 1197 !CDIR NODEP … … 1298 1238 END DO ! jl2 (new ridges) 1299 1239 1300 DO jl2 = ice_cat_bounds(1,1), ice_cat_bounds(1,2)1240 DO jl2 = 1, jpl 1301 1241 1302 1242 !CDIR NODEP … … 1361 1301 CALL wrk_dealloc( jpi, jpj, vrdg1, vrdg2, vsw , srdg1, srdg2, smsw ) 1362 1302 CALL wrk_dealloc( jpi, jpj, afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 ) 1363 CALL wrk_dealloc( jpi, jpj, jpl, aicen_init, vicen_init, vsn on_init, esnon_init, smv_i_init, oa_i_init )1364 CALL wrk_dealloc( jpi, jpj, jkmax, eirft, erdg1, erdg2, ersw )1365 CALL wrk_dealloc( jpi, jpj, jkmax, jpl, eicen_init )1303 CALL wrk_dealloc( jpi, jpj, jpl, aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init ) 1304 CALL wrk_dealloc( jpi, jpj, nlay_i+1, eirft, erdg1, erdg2, ersw ) 1305 CALL wrk_dealloc( jpi, jpj, nlay_i+1, jpl, eicen_init ) 1366 1306 ! 1367 1307 END SUBROUTINE lim_itd_me_ridgeshift … … 1404 1344 !!------------------------------------------------------------------- 1405 1345 INTEGER :: ios ! Local integer output status for namelist read 1406 NAMELIST/namiceitdme/ ridge_scheme_swi, Cs, Cf, fsnowrdg, fsnowrft,& 1407 Gstar, astar, & 1408 Hstar, raftswi, hparmeter, Craft, ridge_por, & 1409 sal_max_ridge, partfun_swi, transfun_swi, & 1410 brinstren_swi 1346 NAMELIST/namiceitdme/ ridge_scheme_swi, Cs, Cf, fsnowrdg, fsnowrft, & 1347 & Gstar, astar, Hstar, raft_swi, hparmeter, Craft, ridge_por, & 1348 & partfun_swi, brinstren_swi 1411 1349 !!------------------------------------------------------------------- 1412 1350 ! … … 1432 1370 WRITE(numout,*)' Equivalent to G* for an exponential part function astar ', astar 1433 1371 WRITE(numout,*)' Quantity playing a role in max ridged ice thickness Hstar ', Hstar 1434 WRITE(numout,*)' Rafting of ice sheets or not raft swi ', raftswi1372 WRITE(numout,*)' Rafting of ice sheets or not raft_swi ', raft_swi 1435 1373 WRITE(numout,*)' Parmeter thickness (threshold between ridge-raft) hparmeter ', hparmeter 1436 1374 WRITE(numout,*)' Rafting hyperbolic tangent coefficient Craft ', Craft 1437 1375 WRITE(numout,*)' Initial porosity of ridges ridge_por ', ridge_por 1438 WRITE(numout,*)' Maximum salinity of ridging ice sal_max_ridge ', sal_max_ridge1439 1376 WRITE(numout,*)' Switch for part. function (0) linear (1) exponential partfun_swi ', partfun_swi 1440 WRITE(numout,*)' Switch for tran. function (0) linear (1) exponential transfun_swi ', transfun_swi1441 1377 WRITE(numout,*)' Switch for including brine volume in ice strength comp. brinstren_swi ', brinstren_swi 1442 1378 ENDIF … … 1462 1398 1463 1399 REAL(wp), POINTER, DIMENSION(:,:) :: zmask ! 2D workspace 1464 REAL(wp) :: zmask_glo 1400 REAL(wp) :: zmask_glo, zsal, zvi, zvs, zei, zes 1465 1401 !!gm REAL(wp) :: xtmp ! temporary variable 1466 1402 !!------------------------------------------------------------------- … … 1468 1404 CALL wrk_alloc( jpi, jpj, zmask ) 1469 1405 1406 ! to be sure that at_i is the sum of a_i(jl) 1407 at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 1408 1470 1409 DO jl = 1, jpl 1471 1472 1410 !----------------------------------------------------------------- 1473 1411 ! Count categories to be zapped. 1474 ! Abort model in case of negative area.1475 1412 !----------------------------------------------------------------- 1476 1413 icells = 0 … … 1478 1415 DO jj = 1, jpj 1479 1416 DO ji = 1, jpi 1480 IF( ( a_i(ji,jj,jl) >= -epsi10 .AND. a_i(ji,jj,jl) < 0._wp ) .OR. & 1481 & ( a_i(ji,jj,jl) > 0._wp .AND. a_i(ji,jj,jl) <= epsi10 ) .OR. & 1482 & ( v_i(ji,jj,jl) == 0._wp .AND. a_i(ji,jj,jl) > 0._wp ) .OR. & 1483 & ( v_i(ji,jj,jl) > 0._wp .AND. v_i(ji,jj,jl) <= epsi10 ) ) zmask(ji,jj) = 1._wp 1417 IF( a_i(ji,jj,jl) <= epsi10 .OR. v_i(ji,jj,jl) <= epsi10 .OR. at_i(ji,jj) <= epsi10 ) THEN 1418 zmask(ji,jj) = 1._wp 1419 ENDIF 1484 1420 END DO 1485 1421 END DO … … 1494 1430 DO jj = 1 , jpj 1495 1431 DO ji = 1 , jpi 1496 !!gm xtmp = e_i(ji,jj,jk,jl) / area(ji,jj) * r1_rdtice 1497 !!gm xtmp = xtmp * unit_fac 1498 ! fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp 1432 zei = e_i(ji,jj,jk,jl) 1499 1433 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * ( 1._wp - zmask(ji,jj) ) 1434 t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * ( 1._wp - zmask(ji,jj) ) + rtt * zmask(ji,jj) 1435 ! update exchanges with ocean 1436 hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 1500 1437 END DO 1501 1438 END DO … … 1504 1441 DO jj = 1 , jpj 1505 1442 DO ji = 1 , jpi 1506 1443 1444 zsal = smv_i(ji,jj,jl) 1445 zvi = v_i(ji,jj,jl) 1446 zvs = v_s(ji,jj,jl) 1447 zes = e_s(ji,jj,1,jl) 1507 1448 !----------------------------------------------------------------- 1508 1449 ! Zap snow energy and use ocean heat to melt snow … … 1514 1455 ! fluxes are positive to the ocean 1515 1456 ! here the flux has to be negative for the ocean 1516 !!gm xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) * r1_rdtice1517 ! fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp1518 1519 !!gm xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) * r1_rdtice !RB ???????1520 1521 1457 t_s(ji,jj,1,jl) = rtt * zmask(ji,jj) + t_s(ji,jj,1,jl) * ( 1._wp - zmask(ji,jj) ) 1522 1458 … … 1524 1460 ! zap ice and snow volume, add water and salt to ocean 1525 1461 !----------------------------------------------------------------- 1526 1527 ! xtmp = (rhoi*vicen(i,j,n) + rhos*vsnon(i,j,n)) / dt 1528 ! sfx_res(ji,jj) = sfx_res(ji,jj) + ( sss_m(ji,jj) ) & 1529 ! * rhosn * v_s(ji,jj,jl) * r1_rdtice 1530 ! sfx_res(ji,jj) = sfx_res(ji,jj) + ( sss_m(ji,jj) - sm_i(ji,jj,jl) ) & 1531 ! * rhoic * v_i(ji,jj,jl) * r1_rdtice 1532 ! sfx (i,j) = sfx (i,j) + xtmp 1533 1534 ato_i(ji,jj) = a_i (ji,jj,jl) * zmask(ji,jj) + ato_i(ji,jj) 1462 ato_i(ji,jj) = a_i (ji,jj,jl) * zmask(ji,jj) + ato_i(ji,jj) 1535 1463 a_i (ji,jj,jl) = a_i (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 1536 1464 v_i (ji,jj,jl) = v_i (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) … … 1539 1467 oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 1540 1468 smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 1541 ! 1469 e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * ( 1._wp - zmask(ji,jj) ) 1470 ! additional condition 1471 IF( v_s(ji,jj,jl) <= epsi10 ) THEN 1472 v_s(ji,jj,jl) = 0._wp 1473 e_s(ji,jj,1,jl) = 0._wp 1474 ENDIF 1475 ! update exchanges with ocean 1476 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice 1477 wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice 1478 wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice 1479 hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 1542 1480 END DO 1543 1481 END DO 1544 ! 1545 END DO ! jl 1482 END DO ! jl 1483 1484 ! to be sure that at_i is the sum of a_i(jl) 1485 at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 1546 1486 ! 1547 1487 CALL wrk_dealloc( jpi, jpj, zmask )
Note: See TracChangeset
for help on using the changeset viewer.