Changeset 5134
- Timestamp:
- 2015-03-09T18:27:34+01:00 (8 years ago)
- Location:
- trunk/NEMOGCM/NEMO/LIM_SRC_3
- Files:
-
- 10 edited
-
limitd_me.F90 (modified) (29 diffs)
-
limitd_th.F90 (modified) (17 diffs)
-
limthd.F90 (modified) (6 diffs)
-
limthd_dh.F90 (modified) (9 diffs)
-
limthd_ent.F90 (modified) (1 diff)
-
limthd_lac.F90 (modified) (2 diffs)
-
limtrp.F90 (modified) (12 diffs)
-
limupdate1.F90 (modified) (3 diffs)
-
limupdate2.F90 (modified) (5 diffs)
-
limvar.F90 (modified) (8 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
r5128 r5134 19 19 USE ice ! LIM variables 20 20 USE dom_ice ! LIM domain 21 USE limthd_lac ! LIM22 21 USE limvar ! LIM 23 USE in_out_manager ! I/O manager24 22 USE lbclnk ! lateral boundary condition - MPP exchanges 25 23 USE lib_mpp ! MPP library … … 27 25 USE prtctl ! Print control 28 26 27 USE in_out_manager ! I/O manager 29 28 USE iom ! I/O manager 30 29 USE lib_fortran ! glob_sum … … 158 157 ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons 159 158 !-----------------------------------------------------------------------------! 160 Cp = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0 159 Cp = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0 ! proport const for PE 161 160 ! 162 161 CALL lim_itd_me_ridgeprep ! prepare ridging … … 275 274 ! 3.4 Compute total area of ice plus open water after ridging. 276 275 !-----------------------------------------------------------------------------! 277 278 CALL lim_itd_me_asumr 276 ! This is in general not equal to one because of divergence during transport 277 asum(:,:) = ato_i(:,:) 278 DO jl = 1, jpl 279 asum(:,:) = asum(:,:) + a_i(:,:,jl) 280 END DO 279 281 280 282 ! 3.5 Do we keep on iterating ??? … … 341 343 342 344 ! Check if there is a ridging error 343 DO jj = 1, jpj 344 DO ji = 1, jpi 345 IF( ABS( asum(ji,jj) - kamax) > epsi10 ) THEN ! there is a bug 346 WRITE(numout,*) ' ' 347 WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj) 348 WRITE(numout,*) ' limitd_me ' 349 WRITE(numout,*) ' POINT : ', ji, jj 350 WRITE(numout,*) ' jpl, a_i, athorn ' 351 WRITE(numout,*) 0, ato_i(ji,jj), athorn(ji,jj,0) 352 DO jl = 1, jpl 353 WRITE(numout,*) jl, a_i(ji,jj,jl), athorn(ji,jj,jl) 354 END DO 355 ENDIF 356 END DO 357 END DO 345 IF( lwp ) THEN 346 DO jj = 1, jpj 347 DO ji = 1, jpi 348 IF( ABS( asum(ji,jj) - kamax) > epsi10 ) THEN ! there is a bug 349 WRITE(numout,*) ' ' 350 WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj) 351 WRITE(numout,*) ' limitd_me ' 352 WRITE(numout,*) ' POINT : ', ji, jj 353 WRITE(numout,*) ' jpl, a_i, athorn ' 354 WRITE(numout,*) 0, ato_i(ji,jj), athorn(ji,jj,0) 355 DO jl = 1, jpl 356 WRITE(numout,*) jl, a_i(ji,jj,jl), athorn(ji,jj,jl) 357 END DO 358 ENDIF 359 END DO 360 END DO 361 END IF 358 362 359 363 ! Conservation check … … 364 368 ENDIF 365 369 366 !-----------------------------------------------------------------------------! 367 ! 6) Updating state variables and trend terms (done in limupdate) 368 !-----------------------------------------------------------------------------! 370 ! updates 369 371 CALL lim_var_glo2eqv 370 372 CALL lim_var_zapsmall 371 373 CALL lim_var_agg( 1 ) 372 374 373 374 IF(ln_ctl) THEN ! Control print 375 !-----------------------------------------------------------------------------! 376 ! control prints 377 !-----------------------------------------------------------------------------! 378 IF(ln_ctl) THEN 375 379 CALL prt_ctl_info(' ') 376 380 CALL prt_ctl_info(' - Cell values : ') … … 473 477 ! PE gain from ridging ice 474 478 !---------------------------- 475 strength(ji,jj) = strength(ji,jj) + aridge(ji,jj,jl) /krdg(ji,jj,jl) &479 strength(ji,jj) = strength(ji,jj) + aridge(ji,jj,jl) / krdg(ji,jj,jl) & 476 480 * z1_3 * ( hrmax(ji,jj,jl)**2 + hrmin(ji,jj,jl)**2 + hrmax(ji,jj,jl) * hrmin(ji,jj,jl) ) 477 481 !!(a**3-b**3)/(a-b) = a*a+ab+b*b … … 528 532 DO ji = 2, jpim1 529 533 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > epsi10) THEN ! ice is present 530 zworka(ji,jj) = 4.0 * strength(ji,jj) & 531 & + strength(ji-1,jj) * tmask(ji-1,jj,1) & 532 & + strength(ji+1,jj) * tmask(ji+1,jj,1) & 533 & + strength(ji,jj-1) * tmask(ji,jj-1,1) & 534 & + strength(ji,jj+1) * tmask(ji,jj+1,1) 535 536 zworka(ji,jj) = zworka(ji,jj) / & 537 & ( 4.0 + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) ) 534 zworka(ji,jj) = ( 4.0 * strength(ji,jj) & 535 & + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) & 536 & + strength(ji,jj-1) * tmask(ji,jj-1,1) + strength(ji,jj+1) * tmask(ji,jj+1,1) & 537 & ) / ( 4.0 + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) ) 538 538 ELSE 539 539 zworka(ji,jj) = 0._wp … … 625 625 626 626 ! Compute total area of ice plus open water. 627 CALL lim_itd_me_asumr 628 ! This is in general not equal to one 629 ! because of divergence during transport 627 ! This is in general not equal to one because of divergence during transport 628 asum(:,:) = ato_i(:,:) 629 DO jl = 1, jpl 630 asum(:,:) = asum(:,:) + a_i(:,:,jl) 631 END DO 630 632 631 633 ! Compute cumulative thickness distribution function … … 678 680 DO ji = 1, jpi 679 681 IF( Gsum(ji,jj,jl) < rn_gstar) THEN 680 athorn(ji,jj,jl) = Gstari * ( Gsum(ji,jj,jl)-Gsum(ji,jj,jl-1)) * &681 (2.0 - (Gsum(ji,jj,jl-1)+Gsum(ji,jj,jl))*Gstari)682 athorn(ji,jj,jl) = Gstari * ( Gsum(ji,jj,jl) - Gsum(ji,jj,jl-1) ) * & 683 & ( 2.0 - (Gsum(ji,jj,jl-1) + Gsum(ji,jj,jl) ) * Gstari ) 682 684 ELSEIF (Gsum(ji,jj,jl-1) < rn_gstar) THEN 683 athorn(ji,jj,jl) = Gstari * ( rn_gstar-Gsum(ji,jj,jl-1)) * &684 (2.0 - (Gsum(ji,jj,jl-1)+rn_gstar)*Gstari)685 athorn(ji,jj,jl) = Gstari * ( rn_gstar - Gsum(ji,jj,jl-1) ) * & 686 & ( 2.0 - ( Gsum(ji,jj,jl-1) + rn_gstar ) * Gstari ) 685 687 ELSE 686 688 athorn(ji,jj,jl) = 0.0 … … 693 695 ! 694 696 zdummy = 1._wp / ( 1._wp - EXP(-astari) ) ! precompute exponential terms using Gsum as a work array 695 696 697 DO jl = -1, jpl 697 698 Gsum(:,:,jl) = EXP( -Gsum(:,:,jl) * astari ) * zdummy 698 END DO !jl699 END DO 699 700 DO jl = 0, jpl 700 701 athorn(:,:,jl) = Gsum(:,:,jl-1) - Gsum(:,:,jl) 701 702 END DO 702 703 ! 703 ENDIF ! nn_partfun704 ENDIF 704 705 705 706 IF( ln_rafting ) THEN ! Ridging and rafting ice participation functions … … 729 730 IF( ln_rafting ) THEN 730 731 731 IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) > epsi10 ) THEN732 IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) > epsi10 .AND. lwp ) THEN 732 733 DO jl = 1, jpl 733 734 DO jj = 1, jpj … … 793 794 ENDIF 794 795 795 END DO ! ji796 END DO ! jj797 END DO ! jl796 END DO 797 END DO 798 END DO 798 799 799 800 ! Normalization factor : aksum, ensures mass conservation … … 834 835 INTEGER :: icells ! number of cells with aicen > puny 835 836 REAL(wp) :: hL, hR, farea, ztmelts ! left and right limits of integration 836 REAL(wp) :: zsstK ! SST in Kelvin837 837 838 838 INTEGER , POINTER, DIMENSION(:) :: indxi, indxj ! compressed indices … … 910 910 ato_i(ji,jj) = 0._wp 911 911 ENDIF 912 END DO !jj913 END DO !ji912 END DO 913 END DO 914 914 915 915 ! if negative open water area alert it 916 IF( neg_ato_i ) THEN ! there is a bug916 IF( neg_ato_i .AND. lwp ) THEN ! there is a bug 917 917 DO jj = 1, jpj 918 918 DO ji = 1, jpi … … 921 921 WRITE(numout,*) 'Ridging error: ato_i < 0' 922 922 WRITE(numout,*) 'ato_i : ', ato_i(ji,jj) 923 ENDIF ! ato_i < -epsi10923 ENDIF 924 924 END DO 925 925 END DO … … 937 937 smv_i_init(:,:,jl) = smv_i(:,:,jl) 938 938 oa_i_init (:,:,jl) = oa_i (:,:,jl) 939 END DO !jl939 END DO 940 940 941 941 esnwn_init(:,:,:) = e_s(:,:,1,:) … … 968 968 indxi(icells) = ji 969 969 indxj(icells) = jj 970 ENDIF ! test on a_icen_init971 END DO ! ji972 END DO ! jj970 ENDIF 971 END DO 972 END DO 973 973 974 974 large_afrac = .false. … … 1094 1094 dhr2(ji,jj) = hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1) 1095 1095 1096 END DO ! ij1096 END DO 1097 1097 1098 1098 !-------------------------------------------------------------------- … … 1112 1112 ! enthalpy of the trapped seawater (J/m2, >0) 1113 1113 ! clem: if sst>0, then ersw <0 (is that possible?) 1114 zsstK = sst_m(ji,jj) + rt0 1115 ersw(ji,jj,jk) = - rhoic * vsw(ji,jj) * rcp * ( zsstK - rt0 ) * r1_nlay_i 1114 ersw(ji,jj,jk) = - rhoic * vsw(ji,jj) * rcp * sst_m(ji,jj) * r1_nlay_i 1116 1115 1117 1116 ! heat flux to the ocean … … 1121 1120 erdg2(ji,jj,jk) = erdg1(ji,jj,jk) + ersw(ji,jj,jk) 1122 1121 1123 END DO ! ij1124 END DO !jk1122 END DO 1123 END DO 1125 1124 1126 1125 … … 1131 1130 jj = indxj(ij) 1132 1131 eice_init(ji,jj) = eice_init(ji,jj) + erdg2(ji,jj,jk) - erdg1(ji,jj,jk) 1133 END DO ! ij1134 END DO !jk1132 END DO 1133 END DO 1135 1134 ENDIF 1136 1135 1137 IF( large_afrac ) THEN ! there is a bug1136 IF( large_afrac .AND. lwp ) THEN ! there is a bug 1138 1137 DO ij = 1, icells 1139 1138 ji = indxi(ij) … … 1146 1145 END DO 1147 1146 ENDIF 1148 IF( large_afrft ) THEN ! there is a bug1147 IF( large_afrft .AND. lwp ) THEN ! there is a bug 1149 1148 DO ij = 1, icells 1150 1149 ji = indxi(ij) … … 1172 1171 ! Transfer area, volume, and energy accordingly. 1173 1172 1174 IF( hrmin(ji,jj,jl1) >= hi_max(jl2) .OR. & 1175 hrmax(ji,jj,jl1) <= hi_max(jl2-1) ) THEN 1173 IF( hrmin(ji,jj,jl1) >= hi_max(jl2) .OR. hrmax(ji,jj,jl1) <= hi_max(jl2-1) ) THEN 1176 1174 hL = 0._wp 1177 1175 hR = 0._wp … … 1199 1197 ji = indxi(ij) 1200 1198 jj = indxj(ij) 1201 e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + fvol(ji,jj) *erdg2(ji,jj,jk)1199 e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + fvol(ji,jj) * erdg2(ji,jj,jk) 1202 1200 END DO 1203 1201 END DO … … 1213 1211 ! thickness category jl2, transfer area, volume, and energy accordingly. 1214 1212 ! 1215 IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. & 1216 hraft(ji,jj,jl1) > hi_max(jl2-1) ) THEN 1213 IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) > hi_max(jl2-1) ) THEN 1217 1214 a_i (ji,jj ,jl2) = a_i (ji,jj ,jl2) + arft2 (ji,jj) 1218 1215 v_i (ji,jj ,jl2) = v_i (ji,jj ,jl2) + virft (ji,jj) … … 1230 1227 ji = indxi(ij) 1231 1228 jj = indxj(ij) 1232 IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. & 1233 hraft(ji,jj,jl1) > hi_max(jl2-1) ) THEN 1229 IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) > hi_max(jl2-1) ) THEN 1234 1230 e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + eirft(ji,jj,jk) 1235 1231 ENDIF … … 1271 1267 ! 1272 1268 END SUBROUTINE lim_itd_me_ridgeshift 1273 1274 1275 SUBROUTINE lim_itd_me_asumr1276 !!-----------------------------------------------------------------------------1277 !! *** ROUTINE lim_itd_me_asumr ***1278 !!1279 !! ** Purpose : finds total fractional area1280 !!1281 !! ** Method : Find the total area of ice plus open water in each grid cell.1282 !! This is similar to the aggregate_area subroutine except that the1283 !! total area can be greater than 1, so the open water area is1284 !! included in the sum instead of being computed as a residual.1285 !!-----------------------------------------------------------------------------1286 INTEGER :: jl ! dummy loop index1287 !!-----------------------------------------------------------------------------1288 !1289 asum(:,:) = ato_i(:,:) ! open water1290 DO jl = 1, jpl ! ice categories1291 asum(:,:) = asum(:,:) + a_i(:,:,jl)1292 END DO1293 !1294 END SUBROUTINE lim_itd_me_asumr1295 1296 1269 1297 1270 SUBROUTINE lim_itd_me_init … … 1351 1324 SUBROUTINE lim_itd_me_icestrength 1352 1325 END SUBROUTINE lim_itd_me_icestrength 1353 SUBROUTINE lim_itd_me_sort1354 END SUBROUTINE lim_itd_me_sort1355 1326 SUBROUTINE lim_itd_me_init 1356 1327 END SUBROUTINE lim_itd_me_init -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
r5123 r5134 25 25 USE ice ! LIM-3 variables 26 26 USE limvar ! LIM-3 variables 27 USE limcons ! LIM-3 conservation28 27 USE prtctl ! Print control 29 28 USE in_out_manager ! I/O manager … … 38 37 PUBLIC lim_itd_th_rem 39 38 PUBLIC lim_itd_th_reb 40 PUBLIC lim_itd_fitline41 PUBLIC lim_itd_shiftice42 39 43 40 !!---------------------------------------------------------------------- … … 129 126 DO jj = 1, jpj 130 127 DO ji = 1, jpi 131 rswitch = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i(ji,jj,jl) +epsi10 ) ) !0 if no ice and 1 if yes128 rswitch = MAX( 0.0, SIGN( 1.0, a_i(ji,jj,jl) - epsi10 ) ) !0 if no ice and 1 if yes 132 129 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * rswitch 133 rswitch = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i_b(ji,jj,jl) + epsi10) )!0 if no ice and 1 if yes130 rswitch = MAX( 0.0, SIGN( 1.0, a_i_b(ji,jj,jl) - epsi10) ) !0 if no ice and 1 if yes 134 131 zht_i_b(ji,jj,jl) = v_i_b(ji,jj,jl) / MAX( a_i_b(ji,jj,jl), epsi10 ) * rswitch 135 IF( a_i(ji,jj,jl) > epsi10 ) zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl) 132 !clem IF( a_i(ji,jj,jl) > epsi10 ) zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl) 133 zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl) 136 134 END DO 137 135 END DO … … 228 226 229 227 IF( a_i(ji,jj,kubnd) > epsi10 ) THEN 230 zhbnew(ji,jj,kubnd) = 3._wp * ht_i(ji,jj,kubnd) - 2._wp * zhbnew(ji,jj,kubnd-1)228 zhbnew(ji,jj,kubnd) = MAX( hi_max(kubnd-1), 3._wp * ht_i(ji,jj,kubnd) - 2._wp * zhbnew(ji,jj,kubnd-1) ) 231 229 ELSE 232 230 zhbnew(ji,jj,kubnd) = hi_max(kubnd) … … 235 233 ENDIF 236 234 237 IF( zhbnew(ji,jj,kubnd) < hi_max(kubnd-1) ) zhbnew(ji,jj,kubnd) = hi_max(kubnd-1) 238 239 END DO !jj 240 END DO !jj 235 END DO 236 END DO 241 237 242 238 !----------------------------------------------------------------------------------------------- … … 252 248 ij = nind_j(ji) 253 249 254 !ji255 250 IF( a_i(ii,ij,klbnd) > epsi10 ) THEN 256 251 zdh0 = zdhice(ii,ij,klbnd) !decrease of ice thickness in the lower category 257 ! ji, a_i > epsi10258 252 IF( zdh0 < 0.0 ) THEN !remove area from category 1 259 ! ji, a_i > epsi10; zdh0 < 0260 253 zdh0 = MIN( -zdh0, hi_max(klbnd) ) 261 254 … … 276 269 v_i(ii,ij,klbnd) = a_i(ii,ij,klbnd) * ht_i(ii,ij,klbnd) ! clem-useless ? 277 270 ENDIF 278 ! ji, a_i > epsi10 279 280 ELSE ! if ice accretion 281 ! ji, a_i > epsi10; zdh0 > 0 271 272 ELSE ! if ice accretion ! a_i > epsi10; zdh0 > 0 282 273 zhbnew(ii,ij,klbnd-1) = MIN( zdh0, hi_max(klbnd) ) 283 274 ! zhbnew was 0, and is shifted to the right to account for thin ice … … 285 276 ENDIF ! zdh0 286 277 287 ! a_i > epsi10288 278 ENDIF ! a_i > epsi10 289 279 290 END DO ! ji280 END DO 291 281 292 282 !- 7.3 g(h) for each thickness category … … 359 349 ij = nind_j(ji) 360 350 IF ( a_i(ii,ij,1) > epsi10 .AND. ht_i(ii,ij,1) < rn_himin ) THEN 361 a_i (ii,ij,1)= a_i(ii,ij,1) * ht_i(ii,ij,1) / rn_himin351 a_i (ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) / rn_himin 362 352 ht_i(ii,ij,1) = rn_himin 363 353 ENDIF … … 515 505 END DO 516 506 517 !---------------------------------------------------------------------------------------------- 518 ! 2) Check for daice or dvice out of range, allowing for roundoff error 519 !---------------------------------------------------------------------------------------------- 520 ! Note: zdaice < 0 or zdvice < 0 usually happens when category jl 521 ! has a small area, with h(n) very close to a boundary. Then 522 ! the coefficients of g(h) are large, and the computed daice and 523 ! dvice can be in error. If this happens, it is best to transfer 524 ! either the entire category or nothing at all, depending on which 525 ! side of the boundary hice(n) lies. 526 !----------------------------------------------------------------- 527 DO jl = klbnd, kubnd-1 528 529 zdaice_negative = .false. 530 zdvice_negative = .false. 531 zdaice_greater_aicen = .false. 532 zdvice_greater_vicen = .false. 533 534 DO jj = 1, jpj 535 DO ji = 1, jpi 536 537 IF (zdonor(ji,jj,jl) > 0) THEN 538 jl1 = zdonor(ji,jj,jl) 539 540 IF (zdaice(ji,jj,jl) < 0.0) THEN 541 IF (zdaice(ji,jj,jl) > -epsi10) THEN 542 IF ( ( jl1 == jl .AND. ht_i(ji,jj,jl1) > hi_max(jl) ) .OR. & 543 ( jl1 == jl+1 .AND. ht_i(ji,jj,jl1) <= hi_max(jl) ) ) THEN 544 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category 545 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 546 ELSE 547 zdaice(ji,jj,jl) = 0.0 ! shift no ice 548 zdvice(ji,jj,jl) = 0.0 549 ENDIF 550 ELSE 551 zdaice_negative = .true. 552 ENDIF 553 ENDIF 554 555 IF (zdvice(ji,jj,jl) < 0.0) THEN 556 IF (zdvice(ji,jj,jl) > -epsi10 ) THEN 557 IF ( ( jl1 == jl .AND. ht_i(ji,jj,jl1) > hi_max(jl) ) .OR. & 558 ( jl1 == jl+1 .AND. ht_i(ji,jj,jl1) <= hi_max(jl) ) ) THEN 559 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category 560 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 561 ELSE 562 zdaice(ji,jj,jl) = 0.0 ! shift no ice 563 zdvice(ji,jj,jl) = 0.0 564 ENDIF 565 ELSE 566 zdvice_negative = .true. 567 ENDIF 568 ENDIF 569 570 ! If daice is close to aicen, set daice = aicen. 571 IF (zdaice(ji,jj,jl) > a_i(ji,jj,jl1) - epsi10 ) THEN 572 IF (zdaice(ji,jj,jl) < a_i(ji,jj,jl1)+epsi10) THEN 573 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) 574 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 575 ELSE 576 zdaice_greater_aicen = .true. 577 ENDIF 578 ENDIF 579 580 IF (zdvice(ji,jj,jl) > v_i(ji,jj,jl1)-epsi10) THEN 581 IF (zdvice(ji,jj,jl) < v_i(ji,jj,jl1)+epsi10) THEN 582 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) 583 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 584 ELSE 585 zdvice_greater_vicen = .true. 586 ENDIF 587 ENDIF 588 589 ENDIF ! donor > 0 590 END DO 591 END DO 592 593 END DO 594 507 !clem: I think the following is wrong (if enabled, it creates negative concentration/volume around -epsi10) 508 ! !---------------------------------------------------------------------------------------------- 509 ! ! 2) Check for daice or dvice out of range, allowing for roundoff error 510 ! !---------------------------------------------------------------------------------------------- 511 ! ! Note: zdaice < 0 or zdvice < 0 usually happens when category jl 512 ! ! has a small area, with h(n) very close to a boundary. Then 513 ! ! the coefficients of g(h) are large, and the computed daice and 514 ! ! dvice can be in error. If this happens, it is best to transfer 515 ! ! either the entire category or nothing at all, depending on which 516 ! ! side of the boundary hice(n) lies. 517 ! !----------------------------------------------------------------- 518 ! DO jl = klbnd, kubnd-1 519 ! 520 ! zdaice_negative = .false. 521 ! zdvice_negative = .false. 522 ! zdaice_greater_aicen = .false. 523 ! zdvice_greater_vicen = .false. 524 ! 525 ! DO jj = 1, jpj 526 ! DO ji = 1, jpi 527 ! 528 ! IF (zdonor(ji,jj,jl) > 0) THEN 529 ! jl1 = zdonor(ji,jj,jl) 530 ! 531 ! IF (zdaice(ji,jj,jl) < 0.0) THEN 532 ! IF (zdaice(ji,jj,jl) > -epsi10) THEN 533 ! IF ( ( jl1 == jl .AND. ht_i(ji,jj,jl1) > hi_max(jl) ) .OR. & 534 ! ( jl1 == jl+1 .AND. ht_i(ji,jj,jl1) <= hi_max(jl) ) ) THEN 535 ! zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category 536 ! zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 537 ! ELSE 538 ! zdaice(ji,jj,jl) = 0.0 ! shift no ice 539 ! zdvice(ji,jj,jl) = 0.0 540 ! ENDIF 541 ! ELSE 542 ! zdaice_negative = .true. 543 ! ENDIF 544 ! ENDIF 545 ! 546 ! IF (zdvice(ji,jj,jl) < 0.0) THEN 547 ! IF (zdvice(ji,jj,jl) > -epsi10 ) THEN 548 ! IF ( ( jl1 == jl .AND. ht_i(ji,jj,jl1) > hi_max(jl) ) .OR. & 549 ! ( jl1 == jl+1 .AND. ht_i(ji,jj,jl1) <= hi_max(jl) ) ) THEN 550 ! zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category 551 ! zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 552 ! ELSE 553 ! zdaice(ji,jj,jl) = 0.0 ! shift no ice 554 ! zdvice(ji,jj,jl) = 0.0 555 ! ENDIF 556 ! ELSE 557 ! zdvice_negative = .true. 558 ! ENDIF 559 ! ENDIF 560 ! 561 ! ! If daice is close to aicen, set daice = aicen. 562 ! IF (zdaice(ji,jj,jl) > a_i(ji,jj,jl1) - epsi10 ) THEN 563 ! IF (zdaice(ji,jj,jl) < a_i(ji,jj,jl1)+epsi10) THEN 564 ! zdaice(ji,jj,jl) = a_i(ji,jj,jl1) 565 ! zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 566 ! ELSE 567 ! zdaice_greater_aicen = .true. 568 ! ENDIF 569 ! ENDIF 570 ! 571 ! IF (zdvice(ji,jj,jl) > v_i(ji,jj,jl1)-epsi10) THEN 572 ! IF (zdvice(ji,jj,jl) < v_i(ji,jj,jl1)+epsi10) THEN 573 ! zdaice(ji,jj,jl) = a_i(ji,jj,jl1) 574 ! zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 575 ! ELSE 576 ! zdvice_greater_vicen = .true. 577 ! ENDIF 578 ! ENDIF 579 ! 580 ! ENDIF ! donor > 0 581 ! END DO 582 ! END DO 583 ! 584 ! END DO 585 !clem 595 586 !------------------------------------------------------------------------------- 596 587 ! 3) Transfer volume and energy between categories … … 614 605 615 606 jl1 = zdonor(ii,ij,jl) 616 rswitch = MAX( 0. 0 , SIGN( 1.0 , v_i(ii,ij,jl1) - epsi10 ) )617 zworka(ii,ij) = zdvice(ii,ij,jl) / MAX( v_i(ii,ij,jl1), epsi 10 ) * rswitch607 rswitch = MAX( 0._wp , SIGN( 1._wp , v_i(ii,ij,jl1) - epsi20 ) ) 608 zworka(ii,ij) = zdvice(ii,ij,jl) / MAX( v_i(ii,ij,jl1), epsi20 ) * rswitch 618 609 IF( jl1 == jl) THEN ; jl2 = jl1+1 619 610 ELSE ; jl2 = jl … … 710 701 ht_i(ji,jj,jl) = v_i (ji,jj,jl) / a_i(ji,jj,jl) 711 702 t_su(ji,jj,jl) = zaTsfn(ji,jj,jl) / a_i(ji,jj,jl) 712 rswitch = 1.0 - MAX( 0.0, SIGN( 1.0, -v_s(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes713 703 ELSE 714 704 ht_i(ji,jj,jl) = 0._wp … … 765 755 DO jj = 1, jpj 766 756 DO ji = 1, jpi 767 IF( a_i(ji,jj,jl) > epsi10 ) THEN 768 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl) 769 ELSE 770 ht_i(ji,jj,jl) = 0._wp 771 ENDIF 757 rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) ) 758 ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * rswitch 772 759 END DO 773 760 END DO … … 775 762 776 763 !------------------------------------------------------------------------------ 777 ! 2) Make sure thickness of cat klbnd is at least hi_max(klbnd) 778 !------------------------------------------------------------------------------ 779 DO jj = 1, jpj 780 DO ji = 1, jpi 781 IF( a_i(ji,jj,klbnd) > epsi10 ) THEN 782 IF( ht_i(ji,jj,klbnd) <= hi_max(0) .AND. hi_max(0) > 0._wp ) THEN 783 a_i(ji,jj,klbnd) = v_i(ji,jj,klbnd) / hi_max(0) 784 ht_i(ji,jj,klbnd) = hi_max(0) 785 ENDIF 786 ENDIF 787 END DO 788 END DO 789 790 !------------------------------------------------------------------------------ 791 ! 3) If a category thickness is not in bounds, shift the 764 ! 2) If a category thickness is not in bounds, shift the 792 765 ! entire area, volume, and energy to the neighboring category 793 766 !------------------------------------------------------------------------------ … … 818 791 zdonor(ji,jj,jl) = jl 819 792 ! begin TECLIM change 820 !zdaice(ji,jj,jl) = a_i(ji,jj,jl)821 !zdvice(ji,jj,jl) = v_i(ji,jj,jl)822 793 !zdaice(ji,jj,jl) = a_i(ji,jj,jl) * 0.5_wp 823 794 !zdvice(ji,jj,jl) = v_i(ji,jj,jl)-zdaice(ji,jj,jl)*(hi_max(jl)+hi_max(jl-1)) * 0.5_wp 824 795 ! end TECLIM change 825 796 ! clem: how much of a_i you send in cat sup is somewhat arbitrary 826 zdaice(ji,jj,jl) = a_i(ji,jj,jl) * ( ht_i(ji,jj,jl) - hi_max(jl) + epsi 10 ) / ht_i(ji,jj,jl)827 zdvice(ji,jj,jl) = v_i(ji,jj,jl) - ( a_i(ji,jj,jl) - zdaice(ji,jj,jl) ) * ( hi_max(jl) - epsi 10 )797 zdaice(ji,jj,jl) = a_i(ji,jj,jl) * ( ht_i(ji,jj,jl) - hi_max(jl) + epsi20 ) / ht_i(ji,jj,jl) 798 zdvice(ji,jj,jl) = v_i(ji,jj,jl) - ( a_i(ji,jj,jl) - zdaice(ji,jj,jl) ) * ( hi_max(jl) - epsi20 ) 828 799 ENDIF 829 800 END DO … … 839 810 ENDIF 840 811 ! 841 END DO ! jl812 END DO 842 813 843 814 !---------------------------- … … 888 859 889 860 !------------------------------------------------------------------------------ 890 ! 4) Conservation check861 ! 3) Conservation check 891 862 !------------------------------------------------------------------------------ 892 863 -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
r5128 r5134 115 115 DO ji = 1, jpi 116 116 !0 if no ice and 1 if yes 117 rswitch = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) +epsi20 ) )117 rswitch = MAX( 0._wp , SIGN( 1._wp , v_i(ji,jj,jl) - epsi20 ) ) 118 118 !Energy of melting q(S,T) [J.m-3] 119 119 e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * REAL( nlay_i ) … … 125 125 DO ji = 1, jpi 126 126 !0 if no ice and 1 if yes 127 rswitch = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_s(ji,jj,jl) +epsi20 ) )127 rswitch = MAX( 0._wp , SIGN( 1._wp , v_s(ji,jj,jl) - epsi20 ) ) 128 128 !Energy of melting q(S,T) [J.m-3] 129 129 e_s(ji,jj,jk,jl) = rswitch * e_s(ji,jj,jk,jl) / MAX( v_s(ji,jj,jl) , epsi20 ) * REAL( nlay_s ) … … 158 158 DO jj = 1, jpj 159 159 DO ji = 1, jpi 160 rswitch = tmask(ji,jj,1) * ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - at_i(ji,jj) + epsi10 )) ) ! 0 if no ice160 rswitch = tmask(ji,jj,1) * MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 0 if no ice 161 161 ! 162 162 ! ! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget … … 360 360 CALL lim_var_eqv2glo 361 361 362 CALL lim_var_zapsmall 362 363 !-------------------------------------------- 363 364 ! Diagnostic thermodynamic growth rates … … 413 414 414 415 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 415 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)416 416 !------------------------------------------------------------------------------| 417 417 ! 7) Add frazil ice growing in leads. 418 418 !------------------------------------------------------------------------------| 419 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 419 420 CALL lim_thd_lac 420 421 CALL lim_var_glo2eqv ! only for info … … 513 514 a_i_1d(ji) = MAX( 0._wp, a_i_1d(ji) + zda_mel ) 514 515 ! adjust thickness 515 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - a_i_1d(ji) +epsi20 ) )516 rswitch = MAX( 0._wp , SIGN( 1._wp , a_i_1d(ji) - epsi20 ) ) 516 517 ht_i_1d(ji) = rswitch * zv / MAX( a_i_1d(ji), epsi20 ) 517 518 ! retrieve total concentration -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
r5128 r5134 244 244 ! thickness change 245 245 IF( zdh_s_pre(ji) > 0._wp ) THEN 246 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zqprec(ji) +epsi20 ) )246 rswitch = MAX( 0._wp , SIGN( 1._wp , zqprec(ji) - epsi20 ) ) 247 247 zdh_s_mel (ji) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 ) 248 248 zdh_s_mel (ji) = MAX( - zdh_s_pre(ji), zdh_s_mel(ji) ) ! bound melting … … 266 266 ! thickness change 267 267 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) ) 268 rswitch = rswitch * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, - q_s_1d(ji,jk) +epsi20 ) ) )268 rswitch = rswitch * ( MAX( 0._wp, SIGN( 1._wp, q_s_1d(ji,jk) - epsi20 ) ) ) 269 269 zdeltah (ji,jk) = - rswitch * zq_su(ji) / MAX( q_s_1d(ji,jk), epsi20 ) 270 270 zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji) ) ! bound melting … … 321 321 DO jk = 1, nlay_s 322 322 DO ji = kideb,kiut 323 rswitch = MAX( 0._wp , SIGN( 1._wp, - ht_s_1d(ji) +epsi20 ) )324 q_s_1d(ji,jk) = ( 1._wp - rswitch ) / MAX( ht_s_1d(ji), epsi20 ) *&325 & ( ( MAX( 0._wp, dh_s_tot(ji) ) ) * zqprec(ji) + &323 rswitch = MAX( 0._wp , SIGN( 1._wp, ht_s_1d(ji) - epsi20 ) ) 324 q_s_1d(ji,jk) = rswitch / MAX( ht_s_1d(ji), epsi20 ) * & 325 & ( ( MAX( 0._wp, dh_s_tot(ji) ) ) * zqprec(ji) + & 326 326 & ( - MAX( 0._wp, dh_s_tot(ji) ) + ht_s_1d(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) ) 327 327 zq_s(ji) = zq_s(ji) + q_s_1d(ji,jk) … … 453 453 q_i_1d(ji,nlay_i+1) = -zEi * rhoic ! New ice energy of melting (J/m3, >0) 454 454 455 ENDIF ! fc_bo_i456 END DO ! ji457 END DO ! iter455 ENDIF 456 END DO 457 END DO 458 458 459 459 ! Contribution to Energy and Salt Fluxes … … 592 592 zq_rema(ji) = zq_su(ji) + zq_bo(ji) 593 593 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) ) ! =1 if snow 594 rswitch = rswitch * MAX( 0._wp, SIGN( 1._wp, zq_s(ji) - epsi20 ) )595 zdeltah (ji,1) = - rswitch * zq_rema(ji) / MAX( zq_s(ji), epsi20 )594 rswitch = rswitch * MAX( 0._wp, SIGN( 1._wp, q_s_1d(ji,1) - epsi20 ) ) 595 zdeltah (ji,1) = - rswitch * zq_rema(ji) / MAX( q_s_1d(ji,1), epsi20 ) 596 596 zdeltah (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_1d(ji) ) ) ! bound melting 597 597 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,1) … … 599 599 ht_s_1d (ji) = ht_s_1d(ji) + zdeltah(ji,1) 600 600 601 zq_rema(ji) = zq_rema(ji) + zdeltah(ji,1) * zq_s(ji) ! update available heat (J.m-2)601 zq_rema(ji) = zq_rema(ji) + zdeltah(ji,1) * q_s_1d(ji,1) ! update available heat (J.m-2) 602 602 ! heat used to melt snow 603 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * zq_s(ji) * r1_rdtice ! W.m-2 (>0)603 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * q_s_1d(ji,1) * r1_rdtice ! W.m-2 (>0) 604 604 ! Contribution to mass flux 605 605 wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice … … 630 630 631 631 ! entrapment during snow ice formation 632 ! new salinity difference stored (to be used in limthd_ ent.F90)632 ! new salinity difference stored (to be used in limthd_sal.F90) 633 633 IF ( nn_icesal == 2 ) THEN 634 rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi 10 ) )634 rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi20 ) ) 635 635 ! salinity dif due to snow-ice formation 636 dsm_i_si_1d(ji) = ( zs_snic - sm_i_1d(ji) ) * dh_snowice(ji) / MAX( ht_i_1d(ji), epsi 10 ) * rswitch636 dsm_i_si_1d(ji) = ( zs_snic - sm_i_1d(ji) ) * dh_snowice(ji) / MAX( ht_i_1d(ji), epsi20 ) * rswitch 637 637 ! salinity dif due to bottom growth 638 638 IF ( zf_tt(ji) < 0._wp ) THEN 639 dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_1d(ji) ) * dh_i_bott(ji) / MAX( ht_i_1d(ji), epsi 10 ) * rswitch639 dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_1d(ji) ) * dh_i_bott(ji) / MAX( ht_i_1d(ji), epsi20 ) * rswitch 640 640 ENDIF 641 641 ENDIF … … 663 663 h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji) 664 664 665 ! Total ablation (to debug) 666 IF( ht_i_1d(ji) <= 0._wp ) a_i_1d(ji) = 0._wp 667 668 END DO !ji 665 END DO 669 666 670 667 ! … … 676 673 rswitch = 1.0 - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) ) 677 674 t_su_1d(ji) = rswitch * t_su_1d(ji) + ( 1.0 - rswitch ) * rt0 678 END DO ! ji675 END DO 679 676 680 677 DO jk = 1, nlay_s 681 678 DO ji = kideb,kiut 682 679 ! mask enthalpy 683 rswitch = MAX( 0._wp , SIGN( 1._wp, - ht_s_1d(ji) ) )684 q_s_1d(ji,jk) = ( 1.0 - rswitch )* q_s_1d(ji,jk)680 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp, - ht_s_1d(ji) ) ) 681 q_s_1d(ji,jk) = rswitch * q_s_1d(ji,jk) 685 682 ! recalculate t_s_1d from q_s_1d 686 t_s_1d(ji,jk) = rt0 + ( 1._wp - rswitch ) * ( - q_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic ) 687 END DO 688 END DO 683 t_s_1d(ji,jk) = rt0 + rswitch * ( - q_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic ) 684 END DO 685 END DO 686 687 ! --- ensure that a_i = 0 where ht_i = 0 --- 688 WHERE( ht_i_1d == 0._wp ) a_i_1d = 0._wp 689 689 690 690 CALL wrk_dealloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_rema ) -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90
r5123 r5134 132 132 DO jk1 = 1, nlay_i 133 133 DO ji = kideb, kiut 134 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zhnew(ji) +epsi20 ) )134 rswitch = MAX( 0._wp , SIGN( 1._wp , zhnew(ji) - epsi20 ) ) 135 135 qnew(ji,jk1) = rswitch * ( zqh_cum1(ji,jk1) - zqh_cum1(ji,jk1-1) ) / MAX( zhnew(ji), epsi20 ) 136 136 ENDDO -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
r5128 r5134 130 130 DO ji = 1, jpi 131 131 !Energy of melting q(S,T) [J.m-3] 132 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) +epsi20 ) ) !0 if no ice132 rswitch = MAX( 0._wp , SIGN( 1._wp , v_i(ji,jj,jl) - epsi20 ) ) !0 if no ice 133 133 e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl), epsi20 ) * REAL( nlay_i, wp ) 134 134 END DO … … 482 482 DO jl = 1, jpl 483 483 DO ji = 1, nbpac 484 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_1d(ji,jl) +epsi20 ) ) ! 0 if no ice and 1 if yes484 rswitch = MAX( 0._wp , SIGN( 1._wp , za_i_1d(ji,jl) - epsi20 ) ) ! 0 if no ice and 1 if yes 485 485 zoa_i_1d(ji,jl) = za_b(ji,jl) * zoa_i_1d(ji,jl) / MAX( za_i_1d(ji,jl) , epsi20 ) * rswitch 486 486 END DO -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
r5128 r5134 68 68 CHARACTER(len=80) :: cltmp 69 69 ! 70 REAL(wp), POINTER, DIMENSION(:,:) :: zsm , zs0at71 REAL(wp), POINTER, DIMENSION(:,:,:) :: z s0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi, zzs0e72 REAL(wp), POINTER, DIMENSION(:,:,:) :: z s0ow73 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: z s0e70 REAL(wp), POINTER, DIMENSION(:,:) :: zsm 71 REAL(wp), POINTER, DIMENSION(:,:,:) :: z0ice, z0snw, z0ai, z0es , z0smi , z0oi 72 REAL(wp), POINTER, DIMENSION(:,:,:) :: z0opw 73 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: z0ei 74 74 REAL(wp), POINTER, DIMENSION(:,:,:) :: zviold, zvsold, zsmvold ! old ice volume... 75 75 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhimax ! old ice thickness … … 80 80 IF( nn_timing == 1 ) CALL timing_start('limtrp') 81 81 82 CALL wrk_alloc( jpi,jpj, zsm, z s0at, zatold, zeiold, zesold )83 CALL wrk_alloc( jpi,jpj,jpl, z s0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi, zzs0e)84 CALL wrk_alloc( jpi,jpj,1, z s0ow )85 CALL wrk_alloc( jpi,jpj,nlay_i+1,jpl, z s0e)82 CALL wrk_alloc( jpi,jpj, zsm, zatold, zeiold, zesold ) 83 CALL wrk_alloc( jpi,jpj,jpl, z0ice, z0snw, z0ai, z0es , z0smi , z0oi ) 84 CALL wrk_alloc( jpi,jpj,1, z0opw ) 85 CALL wrk_alloc( jpi,jpj,nlay_i+1,jpl, z0ei ) 86 86 CALL wrk_alloc( jpi,jpj,jpl, zhimax, zviold, zvsold, zsmvold ) 87 87 … … 118 118 ! in case advection creates ice too thick. 119 119 !--------------------------------------------------------------------- 120 zhimax(:,:,:) = ht_i(:,:,:) 120 zhimax(:,:,:) = ht_i(:,:,:) + ht_s(:,:,:) 121 121 DO jl = 1, jpl 122 122 DO jj = 2, jpjm1 123 123 DO ji = 2, jpim1 124 zhimax(ji,jj,jl) = MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) )124 zhimax(ji,jj,jl) = MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) + ht_s(ji-1:ji+1,jj-1:jj+1,jl) ) 125 125 END DO 126 126 END DO … … 155 155 ! transported fields 156 156 !------------------------- 157 z s0ow(:,:,1) = ato_i(:,:) * e12t(:,:)! Open water area158 DO jl = 1, jpl 159 z s0sn (:,:,jl)= v_s (:,:,jl) * e12t(:,:) ! Snow volume160 z s0ice(:,:,jl) = v_i (:,:,jl) * e12t(:,:) ! Ice volume161 z s0a (:,:,jl)= a_i (:,:,jl) * e12t(:,:) ! Ice area162 z s0sm (:,:,jl)= smv_i(:,:,jl) * e12t(:,:) ! Salt content163 z s0oi (:,:,jl) = oa_i (:,:,jl) * e12t(:,:) ! Age content164 z s0c0(:,:,jl) = e_s (:,:,1,jl) * e12t(:,:) ! Snow heat content157 z0opw(:,:,1) = ato_i(:,:) * e12t(:,:) ! Open water area 158 DO jl = 1, jpl 159 z0snw (:,:,jl) = v_s (:,:,jl) * e12t(:,:) ! Snow volume 160 z0ice(:,:,jl) = v_i (:,:,jl) * e12t(:,:) ! Ice volume 161 z0ai (:,:,jl) = a_i (:,:,jl) * e12t(:,:) ! Ice area 162 z0smi (:,:,jl) = smv_i(:,:,jl) * e12t(:,:) ! Salt content 163 z0oi (:,:,jl) = oa_i (:,:,jl) * e12t(:,:) ! Age content 164 z0es (:,:,jl) = e_s (:,:,1,jl) * e12t(:,:) ! Snow heat content 165 165 DO jk = 1, nlay_i 166 z s0e(:,:,jk,jl) = e_i (:,:,jk,jl) * e12t(:,:) ! Ice heat content166 z0ei (:,:,jk,jl) = e_i (:,:,jk,jl) * e12t(:,:) ! Ice heat content 167 167 END DO 168 168 END DO … … 171 171 IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN !== odd ice time step: adv_x then adv_y ==! 172 172 DO jt = 1, initad 173 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z s0ow (:,:,1), sxopw(:,:), & !--- ice open water area173 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0opw (:,:,1), sxopw(:,:), & !--- ice open water area 174 174 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 175 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z s0ow (:,:,1), sxopw(:,:), &175 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0opw (:,:,1), sxopw(:,:), & 176 176 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 177 177 DO jl = 1, jpl 178 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z s0ice(:,:,jl), sxice(:,:,jl), & !--- ice volume ---178 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl), & !--- ice volume --- 179 179 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 180 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z s0ice(:,:,jl), sxice(:,:,jl), &180 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl), & 181 181 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 182 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z s0sn(:,:,jl), sxsn (:,:,jl), & !--- snow volume ---182 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl), & !--- snow volume --- 183 183 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 184 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z s0sn(:,:,jl), sxsn (:,:,jl), &184 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl), & 185 185 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 186 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z s0sm(:,:,jl), sxsal(:,:,jl), & !--- ice salinity ---186 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl), & !--- ice salinity --- 187 187 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 188 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z s0sm(:,:,jl), sxsal(:,:,jl), &188 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl), & 189 189 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 190 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z s0oi(:,:,jl), sxage(:,:,jl), & !--- ice age ---190 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0oi (:,:,jl), sxage(:,:,jl), & !--- ice age --- 191 191 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 192 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z s0oi(:,:,jl), sxage(:,:,jl), &192 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0oi (:,:,jl), sxage(:,:,jl), & 193 193 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 194 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z s0a(:,:,jl), sxa (:,:,jl), & !--- ice concentrations ---194 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0ai (:,:,jl), sxa (:,:,jl), & !--- ice concentrations --- 195 195 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 196 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z s0a(:,:,jl), sxa (:,:,jl), &196 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ai (:,:,jl), sxa (:,:,jl), & 197 197 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 198 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z s0c0(:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents ---198 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0es (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents --- 199 199 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 200 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z s0c0(:,:,jl), sxc0 (:,:,jl), &200 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0es (:,:,jl), sxc0 (:,:,jl), & 201 201 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 202 202 DO jk = 1, nlay_i !--- ice heat contents --- 203 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z s0e(:,:,jk,jl), sxe (:,:,jk,jl), &203 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl), & 204 204 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), & 205 205 & syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 206 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z s0e(:,:,jk,jl), sxe (:,:,jk,jl), &206 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl), & 207 207 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), & 208 208 & syye(:,:,jk,jl), sxye(:,:,jk,jl) ) … … 212 212 ELSE 213 213 DO jt = 1, initad 214 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z s0ow (:,:,1), sxopw(:,:), & !--- ice open water area214 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0opw (:,:,1), sxopw(:,:), & !--- ice open water area 215 215 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 216 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z s0ow (:,:,1), sxopw(:,:), &216 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0opw (:,:,1), sxopw(:,:), & 217 217 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 218 218 DO jl = 1, jpl 219 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z s0ice(:,:,jl), sxice(:,:,jl), & !--- ice volume ---219 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl), & !--- ice volume --- 220 220 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 221 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z s0ice(:,:,jl), sxice(:,:,jl), &221 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl), & 222 222 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 223 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z s0sn(:,:,jl), sxsn (:,:,jl), & !--- snow volume ---223 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl), & !--- snow volume --- 224 224 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 225 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z s0sn(:,:,jl), sxsn (:,:,jl), &225 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl), & 226 226 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 227 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z s0sm(:,:,jl), sxsal(:,:,jl), & !--- ice salinity ---227 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl), & !--- ice salinity --- 228 228 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 229 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z s0sm(:,:,jl), sxsal(:,:,jl), &229 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl), & 230 230 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 231 231 232 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z s0oi(:,:,jl), sxage(:,:,jl), & !--- ice age ---232 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0oi (:,:,jl), sxage(:,:,jl), & !--- ice age --- 233 233 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 234 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z s0oi(:,:,jl), sxage(:,:,jl), &234 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0oi (:,:,jl), sxage(:,:,jl), & 235 235 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 236 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z s0a(:,:,jl), sxa (:,:,jl), & !--- ice concentrations ---236 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0ai (:,:,jl), sxa (:,:,jl), & !--- ice concentrations --- 237 237 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 238 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z s0a(:,:,jl), sxa (:,:,jl), &238 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0ai (:,:,jl), sxa (:,:,jl), & 239 239 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 240 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z s0c0(:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents ---240 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0es (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents --- 241 241 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 242 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z s0c0(:,:,jl), sxc0 (:,:,jl), &242 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0es (:,:,jl), sxc0 (:,:,jl), & 243 243 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 244 244 DO jk = 1, nlay_i !--- ice heat contents --- 245 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z s0e(:,:,jk,jl), sxe (:,:,jk,jl), &245 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl), & 246 246 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), & 247 247 & syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 248 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z s0e(:,:,jk,jl), sxe (:,:,jk,jl), &248 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl), & 249 249 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), & 250 250 & syye(:,:,jk,jl), sxye(:,:,jk,jl) ) … … 257 257 ! Recover the properties from their contents 258 258 !------------------------------------------- 259 ato_i(:,:) = z s0ow(:,:,1) * r1_e12t(:,:)260 DO jl = 1, jpl 261 v_i (:,:,jl) = z s0ice(:,:,jl) * r1_e12t(:,:)262 v_s (:,:,jl) = z s0sn(:,:,jl) * r1_e12t(:,:)263 smv_i(:,:,jl) = z s0sm(:,:,jl) * r1_e12t(:,:)264 oa_i (:,:,jl) = z s0oi (:,:,jl) * r1_e12t(:,:)265 a_i (:,:,jl) = z s0a(:,:,jl) * r1_e12t(:,:)266 e_s (:,:,1,jl) = z s0c0(:,:,jl) * r1_e12t(:,:)259 ato_i(:,:) = z0opw(:,:,1) * r1_e12t(:,:) 260 DO jl = 1, jpl 261 v_i (:,:,jl) = z0ice(:,:,jl) * r1_e12t(:,:) 262 v_s (:,:,jl) = z0snw(:,:,jl) * r1_e12t(:,:) 263 smv_i(:,:,jl) = z0smi(:,:,jl) * r1_e12t(:,:) 264 oa_i (:,:,jl) = z0oi (:,:,jl) * r1_e12t(:,:) 265 a_i (:,:,jl) = z0ai (:,:,jl) * r1_e12t(:,:) 266 e_s (:,:,1,jl) = z0es (:,:,jl) * r1_e12t(:,:) 267 267 DO jk = 1, nlay_i 268 e_i (:,:,jk,jl) = zs0e(:,:,jk,jl) * r1_e12t(:,:)268 e_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e12t(:,:) 269 269 END DO 270 270 END DO … … 293 293 END DO 294 294 ! 295 CALL lim_hdf( ato_i (:,:) ) ! Diffusion295 CALL lim_hdf( ato_i (:,:) ) 296 296 297 297 !------------------------------------ … … 356 356 357 357 IF ( v_i(ji,jj,jl) > 0._wp ) THEN 358 358 359 zvi = v_i (ji,jj,jl) 359 360 zvs = v_s (ji,jj,jl) … … 361 362 zes = e_s (ji,jj,1,jl) 362 363 zei = SUM( e_i(ji,jj,1:nlay_i,jl) ) 363 zdv = v_i(ji,jj,jl) - zviold(ji,jj,jl) 364 365 rswitch = 1._wp 366 IF ( ( zdv > 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) .AND. zatold(ji,jj) < 0.80 ) .OR. & 367 & ( zdv < 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) ) ) THEN 368 ht_i(ji,jj,jl) = MIN( zhimax(ji,jj,jl), hi_max(jl) ) 369 rswitch = MAX( 0._wp, SIGN( 1._wp, ht_i(ji,jj,jl) - epsi20 ) ) 370 a_i(ji,jj,jl) = rswitch * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi20 ) 371 ELSE 372 ht_i(ji,jj,jl) = MAX( MIN( ht_i(ji,jj,jl), hi_max(jl) ), hi_max(jl-1) ) 373 rswitch = MAX( 0._wp, SIGN( 1._wp, ht_i(ji,jj,jl) - epsi20 ) ) 374 a_i(ji,jj,jl) = rswitch * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi20 ) 364 365 zdv = v_i(ji,jj,jl) + v_s(ji,jj,jl) - zviold(ji,jj,jl) - zvsold(ji,jj,jl) 366 367 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. & 368 & ( zdv <= 0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) ) ) THEN 369 370 rswitch = MAX( 0._wp, SIGN( 1._wp, zhimax(ji,jj,jl) - epsi20 ) ) 371 a_i(ji,jj,jl) = rswitch * ( v_i(ji,jj,jl) + v_s(ji,jj,jl) ) / MAX( zhimax(ji,jj,jl), epsi20 ) 372 373 ! small correction due to *rswitch for a_i 374 v_i (ji,jj,jl) = rswitch * v_i (ji,jj,jl) 375 v_s (ji,jj,jl) = rswitch * v_s (ji,jj,jl) 376 smv_i(ji,jj,jl) = rswitch * smv_i(ji,jj,jl) 377 e_s(ji,jj,1,jl) = rswitch * e_s(ji,jj,1,jl) 378 e_i(ji,jj,1:nlay_i,jl) = rswitch * e_i(ji,jj,1:nlay_i,jl) 379 380 ! Update mass fluxes 381 wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice 382 wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice 383 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice 384 hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0 385 hfx_res(ji,jj) = hfx_res(ji,jj) + ( SUM( e_i(ji,jj,1:nlay_i,jl) ) - zei ) * r1_rdtice ! W.m-2 <0 386 375 387 ENDIF 376 388 377 ! small correction due to *rswitch for a_i378 v_i (ji,jj,jl) = rswitch * v_i (ji,jj,jl)379 v_s (ji,jj,jl) = rswitch * v_s (ji,jj,jl)380 smv_i(ji,jj,jl) = rswitch * smv_i(ji,jj,jl)381 e_s(ji,jj,1,jl) = rswitch * e_s(ji,jj,1,jl)382 e_i(ji,jj,1:nlay_i,jl) = rswitch * e_i(ji,jj,1:nlay_i,jl)383 384 ! Update mass fluxes385 wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice386 wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice387 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice388 hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0389 hfx_res(ji,jj) = hfx_res(ji,jj) + ( SUM( e_i(ji,jj,1:nlay_i,jl) ) - zei ) * r1_rdtice ! W.m-2 <0390 389 ENDIF 391 390 … … 423 422 vt_s (:,:) = 0._wp 424 423 at_i (:,:) = 0._wp 425 !426 424 DO jl = 1, jpl 427 425 DO jj = 1, jpj 428 426 DO ji = 1, jpi 429 ! 430 vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl) ! ice volume 431 vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl) ! snow volume 432 at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration 433 END DO 434 END DO 435 END DO 436 ! ------------------------------------------------- 437 438 ! open water 427 vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl) 428 vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl) 429 at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) 430 END DO 431 END DO 432 END DO 433 434 ! --- open water = 1 if at_i=0 -------------------------------- 439 435 DO jj = 1, jpj 440 436 DO ji = 1, jpi 441 ! open water = 1 if at_i=0442 437 rswitch = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) ) 443 438 ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * ato_i(ji,jj) … … 457 452 IF( ln_icectl ) CALL lim_prt( kt, iiceprt, jiceprt,-1, ' - ice dyn & trp - ' ) 458 453 ! 459 CALL wrk_dealloc( jpi,jpj, zsm, z s0at, zatold, zeiold, zesold )460 CALL wrk_dealloc( jpi,jpj,jpl, z s0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi, zzs0e)461 CALL wrk_dealloc( jpi,jpj,1, z s0ow )462 CALL wrk_dealloc( jpi,jpj,nlay_i+1,jpl, z s0e)454 CALL wrk_dealloc( jpi,jpj, zsm, zatold, zeiold, zesold ) 455 CALL wrk_dealloc( jpi,jpj,jpl, z0ice, z0snw, z0ai, z0es , z0smi , z0oi ) 456 CALL wrk_dealloc( jpi,jpj,1, z0opw ) 457 CALL wrk_dealloc( jpi,jpj,nlay_i+1,jpl, z0ei ) 463 458 CALL wrk_dealloc( jpi,jpj,jpl, zviold, zvsold, zhimax, zsmvold ) 464 459 ! -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90
r5123 r5134 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 !-----------------72 ! zap small values73 !-----------------74 CALL lim_var_zapsmall75 76 71 CALL lim_var_glo2eqv 77 78 72 !---------------------------------------------------- 79 ! Rebin categories with thickness out of bounds 80 !---------------------------------------------------- 81 IF ( jpl > 1 ) CALL lim_itd_th_reb(1, jpl) 82 73 ! ice concentration should not exceed amax 74 !----------------------------------------------------- 83 75 at_i(:,:) = 0._wp 84 76 DO jl = 1, jpl … … 86 78 END DO 87 79 88 !----------------------------------------------------89 ! ice concentration should not exceed amax90 !-----------------------------------------------------91 80 DO jl = 1, jpl 92 81 DO jj = 1, jpj … … 94 83 IF( at_i(ji,jj) > rn_amax .AND. a_i(ji,jj,jl) > 0._wp ) THEN 95 84 a_i(ji,jj,jl) = a_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax / at_i(ji,jj) ) ) 96 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl)97 85 ENDIF 98 86 END DO 99 87 END DO 100 88 END DO 101 102 at_i(:,:) = 0._wp103 DO jl = 1, jpl104 at_i(:,:) = a_i(:,:,jl) + at_i(:,:)105 END DO106 89 107 ! 108 ! Final thickness distribution rebinning109 ! 90 !---------------------------------------------------- 91 ! Rebin categories with thickness out of bounds 92 !---------------------------------------------------- 110 93 IF ( jpl > 1 ) CALL lim_itd_th_reb(1, jpl) 111 94 -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90
r5128 r5134 56 56 INTEGER, INTENT(in) :: kt ! number of iteration 57 57 INTEGER :: ji, jj, jk, jl ! dummy loop indices 58 REAL(wp) :: z h, zsal58 REAL(wp) :: zsal 59 59 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 60 60 !!------------------------------------------------------------------- … … 69 69 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate2', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 70 70 71 !-----------------72 ! zap small values73 !-----------------74 CALL lim_var_agg( 1 )75 CALL lim_var_zapsmall76 CALL lim_var_glo2eqv77 78 !----------------------------------------------------79 ! Rebin categories with thickness out of bounds80 !----------------------------------------------------81 IF ( jpl > 1 ) CALL lim_itd_th_reb(1, jpl)82 83 71 !---------------------------------------------------------------------- 84 72 ! Constrain the thickness of the smallest category above himin 85 73 !---------------------------------------------------------------------- 74 CALL lim_var_glo2eqv 86 75 DO jj = 1, jpj 87 76 DO ji = 1, jpi 88 77 IF( v_i(ji,jj,1) > 0._wp .AND. ht_i(ji,jj,1) < rn_himin ) THEN 89 zh = rn_himin / ht_i(ji,jj,1) 90 ht_s(ji,jj,1) = ht_s(ji,jj,1) * zh 91 ht_i(ji,jj,1) = ht_i(ji,jj,1) * zh 92 a_i (ji,jj,1) = a_i(ji,jj,1) / zh 78 a_i (ji,jj,1) = a_i(ji,jj,1) * ht_i(ji,jj,1) / rn_himin 93 79 ENDIF 94 80 END DO … … 108 94 IF( at_i(ji,jj) > rn_amax .AND. a_i(ji,jj,jl) > 0._wp ) THEN 109 95 a_i(ji,jj,jl) = a_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax / at_i(ji,jj) ) ) 110 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl)111 96 ENDIF 112 97 END DO … … 114 99 END DO 115 100 116 at_i(:,:) = 0.0 117 DO jl = 1, jpl 118 at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 119 END DO 120 121 ! -------------------------------------- 122 ! Final thickness distribution rebinning 123 ! -------------------------------------- 101 !---------------------------------------------------- 102 ! Rebin categories with thickness out of bounds 103 !---------------------------------------------------- 124 104 IF ( jpl > 1 ) CALL lim_itd_th_reb( 1, jpl ) 125 105 … … 140 120 ! salinity stays in bounds 141 121 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) ) 142 smv_i(ji,jj,jl) = rswitch * MAX( MIN( rn_simax * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), rn_simin * v_i(ji,jj,jl) ) !+ rn_simin * ( 1._wp - rswitch ) * v_i(ji,jj,jl)122 smv_i(ji,jj,jl) = rswitch * MAX( MIN( rn_simax * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), rn_simin * v_i(ji,jj,jl) ) 143 123 ! associated salt flux 144 124 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice 145 END DO ! ji146 END DO ! jj147 END DO !jl125 END DO 126 END DO 127 END DO 148 128 ENDIF 149 129 -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
r5123 r5134 161 161 DO jj = 1, jpj 162 162 DO ji = 1, jpi 163 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) +epsi10 ) ) !0 if no ice and 1 if yes163 rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) ) !0 if no ice and 1 if yes 164 164 ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * rswitch 165 165 ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * rswitch … … 173 173 DO jj = 1, jpj 174 174 DO ji = 1, jpi 175 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) +epsi10 ) ) !0 if no ice and 1 if yes175 rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) !0 if no ice and 1 if yes 176 176 sm_i(ji,jj,jl) = smv_i(ji,jj,jl) / MAX( v_i(ji,jj,jl) , epsi10 ) * rswitch 177 177 END DO … … 190 190 DO ji = 1, jpi 191 191 ! ! Energy of melting q(S,T) [J.m-3] 192 rswitch = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) +epsi20 ) ) ! rswitch = 0 if no ice and 1 if yes192 rswitch = MAX( 0.0 , SIGN( 1.0 , v_i(ji,jj,jl) - epsi20 ) ) ! rswitch = 0 if no ice and 1 if yes 193 193 zq_i = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * REAL(nlay_i,wp) 194 194 ztmelts = -tmut * s_i(ji,jj,jk,jl) + rt0 ! Ice layer melt temperature … … 215 215 DO ji = 1, jpi 216 216 !Energy of melting q(S,T) [J.m-3] 217 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - v_s(ji,jj,jl) +epsi20 ) ) ! rswitch = 0 if no ice and 1 if yes217 rswitch = MAX( 0._wp , SIGN( 1._wp , v_s(ji,jj,jl) - epsi20 ) ) ! rswitch = 0 if no ice and 1 if yes 218 218 zq_s = rswitch * e_s(ji,jj,jk,jl) / MAX( v_s(ji,jj,jl) , epsi20 ) * REAL(nlay_s,wp) 219 219 ! … … 233 233 DO jj = 1, jpj 234 234 DO ji = 1, jpi 235 rswitch = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) ))235 rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) ) 236 236 tm_i(ji,jj) = tm_i(ji,jj) + rswitch * t_i(ji,jj,jk,jl) * v_i(ji,jj,jl) & 237 237 & / ( REAL(nlay_i,wp) * MAX( vt_i(ji,jj) , epsi10 ) ) … … 339 339 ! ! weighting the profile 340 340 s_i(ji,jj,jk,jl) = zalpha(ji,jj,jl) * zs_zero + ( 1._wp - zalpha(ji,jj,jl) ) * sm_i(ji,jj,jl) 341 END DO ! ji342 END DO ! jj343 END DO ! jk344 END DO ! jl341 END DO 342 END DO 343 END DO 344 END DO 345 345 ! 346 346 ENDIF ! nn_icesal … … 384 384 DO jj = 1, jpj 385 385 DO ji = 1, jpi 386 rswitch = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) ))386 rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) ) 387 387 tm_i(ji,jj) = tm_i(ji,jj) + rswitch * t_i(ji,jj,jk,jl) * v_i(ji,jj,jl) & 388 388 & * r1_nlay_i / MAX( vt_i(ji,jj) , epsi10 ) … … 572 572 smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * rswitch 573 573 574 ! ice salinity must stay in bounds575 IF( nn_icesal == 2 ) THEN576 smv_i(ji,jj,jl) = MAX( MIN( rn_simax * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), rn_simin * v_i(ji,jj,jl) )577 ENDIF578 574 ! update exchanges with ocean 579 575 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
Note: See TracChangeset
for help on using the changeset viewer.
