Changeset 5038 for branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
- Timestamp:
- 2015-01-20T15:26:13+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
r4346 r5038 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 … … 42 42 PUBLIC lim_itd_me_zapsmall 43 43 PUBLIC lim_itd_me_alloc ! called by iceini.F90 44 45 REAL(wp) :: epsi20 = 1.e-20_wp ! constant values46 REAL(wp) :: epsi10 = 1.e-10_wp ! constant values47 REAL(wp) :: epsi06 = 1.e-06_wp ! constant values48 44 49 45 !----------------------------------------------------------------------- … … 143 139 REAL(wp), POINTER, DIMENSION(:,:) :: esnow_mlt ! energy needed to melt snow in ocean (J m-2) 144 140 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... 141 ! 142 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 149 143 !!----------------------------------------------------------------------------- 150 144 IF( nn_timing == 1 ) CALL timing_start('limitd_me') 151 145 152 146 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 147 158 148 IF(ln_ctl) THEN … … 162 152 163 153 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(:,:,:) 154 155 ! conservation test 156 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 179 157 180 158 !-----------------------------------------------------------------------------! … … 362 340 ! 5) Heat, salt and freshwater fluxes 363 341 !-----------------------------------------------------------------------------! 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 ocean342 wfx_snw(ji,jj) = wfx_snw(ji,jj) + msnow_mlt(ji,jj) * r1_rdtice ! fresh water source for ocean 343 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 344 367 345 END DO … … 399 377 CALL lim_itd_me_zapsmall 400 378 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 379 415 380 IF(ln_ctl) THEN ! Control print … … 445 410 ENDIF 446 411 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 ! ------------------------------- 412 ! conservation test 413 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 470 414 471 415 ENDIF ! ln_limdyn=.true. 472 416 ! 473 417 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 418 ! 477 419 IF( nn_timing == 1 ) CALL timing_stop('limitd_me') … … 670 612 !!---------------------------------------------------------------------! 671 613 INTEGER :: ji,jj, jl ! dummy loop indices 672 INTEGER :: krdg_index !673 614 REAL(wp) :: Gstari, astari, hi, hrmean, zdummy ! local scalar 674 615 REAL(wp), POINTER, DIMENSION(:,:) :: zworka ! temporary array used here … … 746 687 !----------------------------------------------------------------- 747 688 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 689 IF( partfun_swi == 0 ) THEN !--- Linear formulation (Thorndike et al., 1975) 690 DO jl = 0, jpl 752 691 DO jj = 1, jpj 753 692 DO ji = 1, jpi … … 772 711 Gsum(:,:,jl) = EXP( -Gsum(:,:,jl) * astari ) * zdummy 773 712 END DO !jl 774 DO jl = 0, ice_cat_bounds(1,2)713 DO jl = 0, jpl 775 714 athorn(:,:,jl) = Gsum(:,:,jl-1) - Gsum(:,:,jl) 776 715 END DO 777 716 ! 778 ENDIF ! krdg_index779 780 IF( raft swi == 1 ) THEN ! Ridging and rafting ice participation functions717 ENDIF ! partfun_swi 718 719 IF( raft_swi == 1 ) THEN ! Ridging and rafting ice participation functions 781 720 ! 782 721 DO jl = 1, jpl … … 794 733 END DO ! jl 795 734 796 ELSE ! raft swi = 0735 ELSE ! raft_swi = 0 797 736 ! 798 737 DO jl = 1, jpl … … 802 741 ENDIF 803 742 804 IF ( raft swi == 1 ) THEN743 IF ( raft_swi == 1 ) THEN 805 744 806 745 IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) .GT. epsi10 ) THEN … … 908 847 INTEGER :: ij ! horizontal index, combines i and j loops 909 848 INTEGER :: icells ! number of cells with aicen > puny 910 REAL(wp) :: zindb, zsrdg2 ! local scalar911 849 REAL(wp) :: hL, hR, farea, zdummy, zdummy0, ztmelts ! left and right limits of integration 850 REAL(wp) :: zsstK ! SST in Kelvin 912 851 913 852 INTEGER , POINTER, DIMENSION(:) :: indxi, indxj ! compressed indices … … 917 856 918 857 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 ridging858 REAL(wp), POINTER, DIMENSION(:,:,:) :: vsnwn_init, esnwn_init ! snow volume & energy before ridging 920 859 REAL(wp), POINTER, DIMENSION(:,:,:) :: smv_i_init, oa_i_init ! ice salinity & age before ridging 921 860 … … 952 891 CALL wrk_alloc( jpi, jpj, vrdg1, vrdg2, vsw , srdg1, srdg2, smsw ) 953 892 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 )893 CALL wrk_alloc( jpi, jpj, jpl, aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init ) 894 CALL wrk_alloc( jpi, jpj, nlay_i+1, eirft, erdg1, erdg2, ersw ) 895 CALL wrk_alloc( jpi, jpj, nlay_i+1, jpl, eicen_init ) 957 896 958 897 ! Conservation check … … 1008 947 aicen_init(:,:,jl) = a_i(:,:,jl) 1009 948 vicen_init(:,:,jl) = v_i(:,:,jl) 1010 vsn on_init(:,:,jl) = v_s(:,:,jl)949 vsnwn_init(:,:,jl) = v_s(:,:,jl) 1011 950 ! 1012 951 smv_i_init(:,:,jl) = smv_i(:,:,jl) … … 1014 953 END DO !jl 1015 954 1016 esn on_init(:,:,:) = e_s(:,:,1,:)955 esnwn_init(:,:,:) = e_s(:,:,1,:) 1017 956 1018 957 DO jl = 1, jpl … … 1091 1030 ! / rafting category n1. 1092 1031 !-------------------------------------------------------------------------- 1093 vrdg1(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj) / ( 1._wp + ridge_por )1032 vrdg1(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj) 1094 1033 vrdg2(ji,jj) = vrdg1(ji,jj) * ( 1. + ridge_por ) 1095 1034 vsw (ji,jj) = vrdg1(ji,jj) * ridge_por 1096 1035 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) 1036 vsrdg(ji,jj) = vsnwn_init(ji,jj,jl1) * afrac(ji,jj) 1037 esrdg(ji,jj) = esnwn_init(ji,jj,jl1) * afrac(ji,jj) 1038 srdg1(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) 1039 srdg2(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) !! MV HC 2014 this line seems useless 1101 1040 1102 1041 ! rafting volumes, heat contents ... 1103 1042 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)1043 vsrft(ji,jj) = vsnwn_init(ji,jj,jl1) * afrft(ji,jj) 1044 esrft(ji,jj) = esnwn_init(ji,jj,jl1) * afrft(ji,jj) 1106 1045 smrft(ji,jj) = smv_i_init(ji,jj,jl1) * afrft(ji,jj) 1107 1046 … … 1120 1059 ! Salinity 1121 1060 !------------- 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 1061 smsw(ji,jj) = vsw(ji,jj) * sss_m(ji,jj) ! salt content of seawater frozen in voids !! MV HC2014 1062 srdg2(ji,jj) = srdg1(ji,jj) + smsw(ji,jj) ! salt content of new ridge 1063 1064 !srdg2(ji,jj) = MIN( s_i_max * vrdg2(ji,jj) , zsrdg2 ) ! impose a maximum salinity 1127 1065 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 1066 sfx_dyn(ji,jj) = sfx_dyn(ji,jj) - smsw(ji,jj) * rhoic * r1_rdtice 1067 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 1068 1133 1069 !------------------------------------ … … 1158 1094 & + rhosn*vsrft(ji,jj)*(1.0-fsnowrft) 1159 1095 1160 esnow_mlt(ji,jj) = esnow_mlt(ji,jj) + esrdg(ji,jj)*(1.0-fsnowrdg) & !rafting included 1161 & + esrft(ji,jj)*(1.0-fsnowrft) 1096 ! in 1e-9 Joules (same as e_s) 1097 esnow_mlt(ji,jj) = esnow_mlt(ji,jj) - esrdg(ji,jj)*(1.0-fsnowrdg) & !rafting included 1098 & - esrft(ji,jj)*(1.0-fsnowrft) 1162 1099 1163 1100 !----------------------------------------------------------------- … … 1184 1121 jj = indxj(ij) 1185 1122 ! heat content of ridged ice 1186 erdg1(ji,jj,jk) = eicen_init(ji,jj,jk,jl1) * afrac(ji,jj) / ( 1._wp + ridge_por )1123 erdg1(ji,jj,jk) = eicen_init(ji,jj,jk,jl1) * afrac(ji,jj) 1187 1124 eirft(ji,jj,jk) = eicen_init(ji,jj,jk,jl1) * afrft(ji,jj) 1188 1125 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 1126 1127 1128 ! enthalpy of the trapped seawater (J/m2, >0) 1129 ! clem: if sst>0, then ersw <0 (is that possible?) 1130 zsstK = sst_m(ji,jj) + rt0 1131 ersw(ji,jj,jk) = - rhoic * vsw(ji,jj) * rcp * ( zsstK - rt0 ) / REAL( nlay_i ) 1132 1133 ! heat flux to the ocean 1134 hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ersw(ji,jj,jk) * r1_rdtice ! > 0 [W.m-2] ocean->ice flux 1203 1135 1204 1136 ! 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 ) 1137 ersw(ji,jj,jk) = ersw(ji,jj,jk) / unit_fac 1138 1139 ! Mutliply by ice volume, and divide by number of layers to get heat content in 1.e9 J 1140 ! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean 1141 !! MV HC 2014 1142 ersw (ji,jj,jk) = ersw(ji,jj,jk) * area(ji,jj) 1209 1143 1210 1144 erdg2(ji,jj,jk) = erdg1(ji,jj,jk) + ersw(ji,jj,jk) 1145 1211 1146 END DO ! ij 1212 1147 END DO !jk … … 1253 1188 !------------------------------------------------------------------------------- 1254 1189 ! jl1 looping 1-jpl 1255 DO jl2 = ice_cat_bounds(1,1), ice_cat_bounds(1,2)1190 DO jl2 = 1, jpl 1256 1191 ! over categories to which ridged ice is transferred 1257 1192 !CDIR NODEP … … 1298 1233 END DO ! jl2 (new ridges) 1299 1234 1300 DO jl2 = ice_cat_bounds(1,1), ice_cat_bounds(1,2)1235 DO jl2 = 1, jpl 1301 1236 1302 1237 !CDIR NODEP … … 1361 1296 CALL wrk_dealloc( jpi, jpj, vrdg1, vrdg2, vsw , srdg1, srdg2, smsw ) 1362 1297 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 )1298 CALL wrk_dealloc( jpi, jpj, jpl, aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init ) 1299 CALL wrk_dealloc( jpi, jpj, nlay_i+1, eirft, erdg1, erdg2, ersw ) 1300 CALL wrk_dealloc( jpi, jpj, nlay_i+1, jpl, eicen_init ) 1366 1301 ! 1367 1302 END SUBROUTINE lim_itd_me_ridgeshift … … 1404 1339 !!------------------------------------------------------------------- 1405 1340 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 1341 NAMELIST/namiceitdme/ ridge_scheme_swi, Cs, Cf, fsnowrdg, fsnowrft, & 1342 & Gstar, astar, Hstar, raft_swi, hparmeter, Craft, ridge_por, & 1343 & partfun_swi, brinstren_swi 1411 1344 !!------------------------------------------------------------------- 1412 1345 ! … … 1418 1351 READ ( numnam_ice_cfg, namiceitdme, IOSTAT = ios, ERR = 902 ) 1419 1352 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitdme in configuration namelist', lwp ) 1420 WRITE ( numoni, namiceitdme )1353 IF(lwm) WRITE ( numoni, namiceitdme ) 1421 1354 ! 1422 1355 IF (lwp) THEN ! control print … … 1432 1365 WRITE(numout,*)' Equivalent to G* for an exponential part function astar ', astar 1433 1366 WRITE(numout,*)' Quantity playing a role in max ridged ice thickness Hstar ', Hstar 1434 WRITE(numout,*)' Rafting of ice sheets or not raft swi ', raftswi1367 WRITE(numout,*)' Rafting of ice sheets or not raft_swi ', raft_swi 1435 1368 WRITE(numout,*)' Parmeter thickness (threshold between ridge-raft) hparmeter ', hparmeter 1436 1369 WRITE(numout,*)' Rafting hyperbolic tangent coefficient Craft ', Craft 1437 1370 WRITE(numout,*)' Initial porosity of ridges ridge_por ', ridge_por 1438 WRITE(numout,*)' Maximum salinity of ridging ice sal_max_ridge ', sal_max_ridge1439 1371 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 1372 WRITE(numout,*)' Switch for including brine volume in ice strength comp. brinstren_swi ', brinstren_swi 1442 1373 ENDIF … … 1462 1393 1463 1394 REAL(wp), POINTER, DIMENSION(:,:) :: zmask ! 2D workspace 1464 REAL(wp) :: zmask_glo 1395 REAL(wp) :: zmask_glo, zsal, zvi, zvs, zei, zes 1465 1396 !!gm REAL(wp) :: xtmp ! temporary variable 1466 1397 !!------------------------------------------------------------------- … … 1468 1399 CALL wrk_alloc( jpi, jpj, zmask ) 1469 1400 1401 ! to be sure that at_i is the sum of a_i(jl) 1402 at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 1403 1470 1404 DO jl = 1, jpl 1471 1472 1405 !----------------------------------------------------------------- 1473 1406 ! Count categories to be zapped. 1474 ! Abort model in case of negative area.1475 1407 !----------------------------------------------------------------- 1476 1408 icells = 0 … … 1478 1410 DO jj = 1, jpj 1479 1411 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 1412 IF( a_i(ji,jj,jl) <= epsi10 .OR. v_i(ji,jj,jl) <= epsi10 .OR. at_i(ji,jj) <= epsi10 ) THEN 1413 zmask(ji,jj) = 1._wp 1414 ENDIF 1484 1415 END DO 1485 1416 END DO … … 1494 1425 DO jj = 1 , jpj 1495 1426 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 1427 zei = e_i(ji,jj,jk,jl) 1499 1428 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * ( 1._wp - zmask(ji,jj) ) 1429 t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * ( 1._wp - zmask(ji,jj) ) + rtt * zmask(ji,jj) 1430 ! update exchanges with ocean 1431 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 1432 END DO 1501 1433 END DO … … 1504 1436 DO jj = 1 , jpj 1505 1437 DO ji = 1 , jpi 1506 1438 1439 zsal = smv_i(ji,jj,jl) 1440 zvi = v_i(ji,jj,jl) 1441 zvs = v_s(ji,jj,jl) 1442 zes = e_s(ji,jj,1,jl) 1507 1443 !----------------------------------------------------------------- 1508 1444 ! Zap snow energy and use ocean heat to melt snow … … 1514 1450 ! fluxes are positive to the ocean 1515 1451 ! 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 1452 t_s(ji,jj,1,jl) = rtt * zmask(ji,jj) + t_s(ji,jj,1,jl) * ( 1._wp - zmask(ji,jj) ) 1522 1453 … … 1524 1455 ! zap ice and snow volume, add water and salt to ocean 1525 1456 !----------------------------------------------------------------- 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) 1457 ato_i(ji,jj) = a_i (ji,jj,jl) * zmask(ji,jj) + ato_i(ji,jj) 1535 1458 a_i (ji,jj,jl) = a_i (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 1536 1459 v_i (ji,jj,jl) = v_i (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) … … 1539 1462 oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 1540 1463 smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 1541 ! 1464 e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * ( 1._wp - zmask(ji,jj) ) 1465 ! additional condition 1466 IF( v_s(ji,jj,jl) <= epsi10 ) THEN 1467 v_s(ji,jj,jl) = 0._wp 1468 e_s(ji,jj,1,jl) = 0._wp 1469 ENDIF 1470 ! update exchanges with ocean 1471 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice 1472 wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice 1473 wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice 1474 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 1475 END DO 1543 1476 END DO 1544 ! 1545 END DO ! jl 1477 END DO ! jl 1478 1479 ! to be sure that at_i is the sum of a_i(jl) 1480 at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 1546 1481 ! 1547 1482 CALL wrk_dealloc( jpi, jpj, zmask )
Note: See TracChangeset
for help on using the changeset viewer.