Changeset 834 for trunk/NEMO/LIM_SRC_3/limitd_me.F90
- Timestamp:
- 2008-03-07T18:11:35+01:00 (16 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/LIM_SRC_3/limitd_me.F90 (modified) (67 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/LIM_SRC_3/limitd_me.F90
r825 r834 2 2 3 3 #if defined key_lim3 4 !!---------------------------------------------------------------------- 5 !! 'key_lim3' : LIM3 sea-ice model 6 !!---------------------------------------------------------------------- 7 !! 4 8 !!====================================================================== 5 9 !! *** MODULE limitd_me *** … … 7 11 !! computation of changes in g(h) 8 12 !!====================================================================== 9 13 !! 10 14 !!---------------------------------------------------------------------- 11 15 !! * Modules used … … 17 21 USE ice_oce ! ice variables 18 22 USE thd_ice 19 USE limicepoints20 23 USE limistate 21 24 USE in_out_manager 22 25 USE ice 23 26 USE par_ice 24 USE limthd_lab25 27 USE limthd_lac 26 28 USE limvar … … 74 76 REAL(wp) :: & 75 77 Cp 76 77 78 ! 78 79 !----------------------------------------------------------------------- 79 ! 80 ! Ridging diagnostic arrays for history files 80 81 !----------------------------------------------------------------------- 81 82 ! … … 88 89 89 90 !!---------------------------------------------------------------------- 90 !! LIM-@ 3.0, UCL-ASTR (200 5-2006)91 !! LIM-@ 3.0, UCL-ASTR (2008) 91 92 !! (c) UCL-ASTR and Martin Vancoppenolle 92 93 !!---------------------------------------------------------------------- … … 114 115 !! ** External : 115 116 !! 117 !! ** Steps : 118 !! 1) Thickness categories boundaries, ice / o.w. concentrations 119 !! Ridge preparation 120 !! 2) Dynamical inputs (closing rate, divu_adv, opning) 121 !! 3) Ridging iteration 122 !! 4) Ridging diagnostics 123 !! 5) Heat, salt and freshwater fluxes 124 !! 6) Compute increments of tate variables and come back to old values 125 !! 116 126 !! ** References : There are a lot of references and can be difficult / 117 127 !! boring to read … … 138 148 !! 139 149 !! ** History : 140 !! This routine was freely inpired from CICE 141 !! authors William H. Lipscomb, LANL 142 !! Elizabeth C. Hunke (LANL) 143 !! the routines were available online 144 !! so, big thanks to the authors and to LANL 145 !! The modified features are the ridging/rafting algorithm 150 !! This routine is based on CICE code 151 !! and authors William H. Lipscomb, 152 !! and Elizabeth C. Hunke, LANL 153 !! are gratefully acknowledged 146 154 !! 147 155 !! (02-2006) Martin Vancoppenolle, UCL-ASTR 148 !! Voila, c'est la derniere etape du codage!!! J'en ai149 !! plein le cul mais je m'en sors mieux qu'avant!!!150 !! je ne sais pas quand ca va marcher151 !! O toi, futur utilisateur, si tu pestes sur ce code152 !! Pardonne moi, sache que j'ai bien peste dessus aussi153 156 !! 154 157 !!--------------------------------------------------------------------! 155 158 !! * Arguments 156 159 157 !! * Local variables 158 INTEGER :: ji, & ! spatial dummy loop index 159 jj, & ! spatial dummy loop index 160 jk, & ! vertical layering dummy loop index 161 jl, & ! ice category dummy loop index 162 jm, & ! ice types dummy loop index 163 index, & ! indexes for ice points debugging 164 jbnd1, & 165 jbnd2, & 166 niter, & ! iteration counter 167 nitermax = 20 ! max number of ridging iterations (it used to 168 ! be 20 169 170 REAL(wp) :: & ! constant values 171 zeps = 1.0e-10, & 172 epsi10 = 1.0e-10, & 173 epsi06 = 1.0e-6 174 175 REAL(wp) :: & ! constant values for ice enthalpy 176 zindb 177 178 REAL(wp), DIMENSION(jpi,jpj) :: & 179 closing_net, & ! net rate at which area is removed (1/s) 180 ! (ridging ice area - area of new ridges) / dt 181 divu_adv , & ! divu as implied by transport scheme (1/s) 182 opning , & ! rate of opening due to divergence/shear 183 closing_gross, & ! rate at which area removed, not counting 184 ! area of new ridges 185 msnow_mlt , & ! mass of snow added to ocean (kg m-2) 186 esnow_mlt ! energy needed to melt snow in ocean (J m-2) 187 188 REAL(wp) :: & 189 w1, & ! temporary variable 190 tmpfac, & ! factor by which opening/closing rates are cut 191 dti ! 1 / dt 192 193 LOGICAL :: & 194 iterate_ridging, & ! if true, repeat the ridging 195 asum_error ! flag for asum .ne. 1 196 197 REAL(wp) :: & 198 big = 1.0e8 199 200 REAL (wp), DIMENSION(jpi,jpj) :: & ! 201 vt_i_init, vt_i_final, & ! ice volume summed over categories 202 vt_s_init, vt_s_final, & ! snow volume summed over categories 203 et_i_init, et_i_final, & ! ice energy summed over categories 204 et_s_init, et_s_final ! snow energy summed over categories 205 206 CHARACTER (len = 15) :: fieldid 160 !! * Local variables 161 INTEGER :: ji, & ! spatial dummy loop index 162 jj, & ! spatial dummy loop index 163 jk, & ! vertical layering dummy loop index 164 jl, & ! ice category dummy loop index 165 niter, & ! iteration counter 166 nitermax = 20 ! max number of ridging iterations 167 168 REAL(wp) :: & ! constant values 169 zeps = 1.0e-10, & 170 epsi10 = 1.0e-10, & 171 epsi06 = 1.0e-6 172 173 REAL(wp), DIMENSION(jpi,jpj) :: & 174 closing_net, & ! net rate at which area is removed (1/s) 175 ! (ridging ice area - area of new ridges) / dt 176 divu_adv , & ! divu as implied by transport scheme (1/s) 177 opning , & ! rate of opening due to divergence/shear 178 closing_gross, & ! rate at which area removed, not counting 179 ! area of new ridges 180 msnow_mlt , & ! mass of snow added to ocean (kg m-2) 181 esnow_mlt ! energy needed to melt snow in ocean (J m-2) 182 183 REAL(wp) :: & 184 w1, & ! temporary variable 185 tmpfac, & ! factor by which opening/closing rates are cut 186 dti ! 1 / dt 187 188 LOGICAL :: & 189 iterate_ridging, & ! if true, repeat the ridging 190 asum_error ! flag for asum .ne. 1 191 192 REAL(wp) :: & 193 big = 1.0e8 194 195 REAL (wp), DIMENSION(jpi,jpj) :: & ! 196 vt_i_init, vt_i_final ! ice volume summed over categories 197 198 CHARACTER (len = 15) :: fieldid 207 199 208 200 !!-- End of declarations … … 223 215 ! Set hi_max(ncat) to a big value to ensure that all ridged ice 224 216 ! is thinner than hi_max(ncat). 217 225 218 Cp = 0.5* grav * (rau0-rhoic)*rhoic/rau0 ! proport const for PE 226 ! BACK TO LIMA_MEC 227 ! hi_max(ice_cat_bounds(2,2)) = big 228 229 CALL lim_itd_me_ridgeprep 230 231 ! ! conservation check 232 ! CALL lim_column_sum (jpl, v_i, vt_i_init) 219 CALL lim_itd_me_ridgeprep ! prepare ridging 220 221 ! conservation check 222 IF ( con_i) CALL lim_column_sum (jpl, v_i, vt_i_init) 233 223 234 224 ! Initialize arrays. … … 288 278 END DO 289 279 290 !+++++ [291 WRITE(numout,*) ' 2. Dynamical inputs '292 WRITE(numout,*) ' closing_net ', closing_net(jiindex,jjindex)293 WRITE(numout,*) ' divu_adv ', divu_adv (jiindex,jjindex)294 WRITE(numout,*) ' opning ', opning (jiindex,jjindex)295 !+++++ ]296 297 298 280 !-----------------------------------------------------------------------------! 299 281 ! 3) Ridging iteration … … 304 286 305 287 DO WHILE ( iterate_ridging .AND. niter < nitermax ) 306 307 !+++++ [308 ! WRITE(numout,*)309 WRITE(numout,*) ' 3. Ridging iteration : ', niter310 WRITE(numout,*) ' ---------------------- '311 ! WRITE(numout,*)312 ! WRITE(numout,*) 'a_i ', a_i(jiindex, jjindex,1:jpl)313 ! WRITE(numout,*) 'v_i ', v_i(jiindex, jjindex,1:jpl)314 !+++++ ]315 316 ! 3.1 ice / o.w. concentrationts, participation and transfer functions317 288 318 289 DO jj = 1, jpj … … 342 313 END DO ! jj 343 314 344 !+++++ [345 WRITE(numout,*) ' correction 1'346 WRITE(numout,*) ' closing_gross ', closing_gross(jiindex,jjindex)347 WRITE(numout,*) ' opning ', opning (jiindex,jjindex)348 !+++++ [349 350 315 ! correction to closing rate / opening if excessive ice removal 351 316 !--------------------------------------------------------------- … … 356 321 DO jj = 1, jpj 357 322 DO ji = 1, jpi 358 ! IF (numit.GE.565) WRITE(numout,*) 'Stop Hier : ',ji,jj,jl359 323 IF ( a_i(ji,jj,jl) .GT. epsi11 .AND. athorn(ji,jj,jl) .GT. 0.0 ) THEN 360 324 w1 = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice … … 369 333 END DO !jl 370 334 371 !+++++ [372 WRITE(numout,*) ' correction 2'373 WRITE(numout,*) ' closing_gross ', closing_gross(jiindex,jjindex)374 WRITE(numout,*) ' opning ', opning(jiindex,jjindex)375 !+++++ ]376 377 335 ! 3.3 Redistribute area, volume, and energy. 378 336 !-----------------------------------------------------------------------------! … … 385 343 386 344 CALL lim_itd_me_asumr 387 388 !+++++ [389 WRITE(numout,*)390 WRITE(numout,*) ' back from asumr'391 WRITE(numout,*) ' asum, ato_i, at_i : ', asum(jiindex,jjindex), &392 ato_i(jiindex,jjindex), at_i(jiindex, jjindex)393 WRITE(numout,*) ' a_i : ', a_i(jiindex,jjindex,1:jpl)394 WRITE(numout,*) ' v_i : ', v_i(jiindex,jjindex,1:jpl)395 !+++++ ]396 345 397 346 ! 3.5 Do we keep on iterating ??? … … 457 406 !-----------------------------------------------------------------------------! 458 407 ! fresh water source for ocean 459 fmmec(ji,jj) = fmmec(ji,jj) + msnow_mlt(ji,jj)*dti !in kg / m2 /s !!! very important 460 ! fresh_hist(ji,jj) = fresh_hist(ji,jj) + msnow_mlt(ji,jj)*dti 408 fmmec(ji,jj) = fmmec(ji,jj) + msnow_mlt(ji,jj)*dti 461 409 462 410 ! heat sink for ocean 463 411 fhmec(ji,jj) = fhmec(ji,jj) + esnow_mlt(ji,jj)*dti 464 ! fhnet_hist(ji,jj) = fhnet_hist(ji,jj) + esnow_mlt(ji,jj)*dti465 412 466 413 END DO … … 486 433 487 434 ! Conservation check 488 ! CALL lim_column_sum (jpl, v_i, vt_i_final) 489 ! fieldid = ' v_i : limitd_me ' 490 ! CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid) 435 IF ( con_i ) THEN 436 CALL lim_column_sum (jpl, v_i, vt_i_final) 437 fieldid = ' v_i : limitd_me ' 438 CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid) 439 ENDIF 491 440 492 441 !-----------------------------------------------------------------------------! … … 494 443 !-----------------------------------------------------------------------------! 495 444 496 !+++++ 497 ! BACK TO LIMA_MEC 498 ! hi_max(ice_cat_bounds(2,2)) = 20.0 499 ! only for debugging 500 CALL lim_var_glo2eqv 501 502 WRITE(numout,*) ' End of lim_itd_me ****** ' 503 WRITE(numout,*) ' u_ice : ', u_ice(jiindex,jjindex) 504 WRITE(numout,*) ' v_ice : ', v_ice(jiindex,jjindex) 505 506 DO jl = 1, jpl 507 WRITE(numout,*) '* - category number ', jl 508 WRITE(numout,*) ' a_i : ', a_i(jiindex,jjindex,jl) 509 WRITE(numout,*) ' ht_i : ', ht_i(jiindex,jjindex,jl) 510 WRITE(numout,*) ' v_i : ', v_i(jiindex,jjindex,jl) 511 WRITE(numout,*) ' v_s : ', v_s(jiindex,jjindex,jl) 512 WRITE(numout,*) ' e_s : ', e_s(jiindex,jjindex,1,jl)/1.0e9 513 WRITE(numout,*) ' e_i1 : ', e_i(jiindex,jjindex,1,jl)/1.0e9 514 WRITE(numout,*) ' e_i2 : ', e_i(jiindex,jjindex,2,jl)/1.0e9 515 WRITE(numout,*) ' smv_i : ', smv_i(jiindex,jjindex,jl) 516 WRITE(numout,*) ' sm_i : ', sm_i(jiindex,jjindex,jl) 517 WRITE(numout,*) ' oa_i : ', oa_i (jiindex,jjindex,jl) 518 WRITE(numout,*) ' o_i : ', o_i (jiindex,jjindex,jl) 519 WRITE(numout,*) ' t_snow : ', t_s(jiindex,jjindex,1,jl) 520 WRITE(numout,*) ' t_su : ', t_su(jiindex,jjindex,jl) 521 WRITE(numout,*) ' t_i1 : ', t_i(jiindex,jjindex,1,jl) 522 WRITE(numout,*) ' t_i2 : ', t_i(jiindex,jjindex,2,jl) 523 WRITE(numout,*) 524 END DO 525 526 ! DO jl = 1, jpl 527 ! DO jj = 1, jpj 528 ! DO ji = 1, jpi 529 ! IF (old_a_i(ji,jj,jl).eq.0.00) THEN 530 ! o_i(ji,jj,jl) = MAX( MIN( rdt_ice*float(numit)/86400.0, o_i(ji,jj,jl) ), 0.0 ) 531 ! ENDIF 532 ! END DO 533 ! END DO 534 ! END DO 535 !+++++ 536 537 ! ! Zero out categories with very small areas 445 CALL lim_var_glo2eqv 538 446 CALL lim_itd_me_zapsmall 539 447 … … 615 523 END DO 616 524 END DO 617 618 619 525 620 526 END SUBROUTINE lim_itd_me … … 624 530 SUBROUTINE lim_itd_me_icestrength (kstrngth) ! (subroutine 2/6) 625 531 626 !!---------------------------------------------------------------------- -------532 !!---------------------------------------------------------------------- 627 533 !! *** ROUTINE lim_itd_me_icestrength *** 628 534 !! ** Purpose : … … 634 540 !! dissipated per unit area removed from the ice pack under compression, 635 541 !! and assumed proportional to the change in potential energy caused 636 !! by ridging. 542 !! by ridging. Note that only Hibler's formulation is stable and that 543 !! ice strength has to be smoothed 637 544 !! 638 545 !! ** Inputs / Ouputs : kstrngth (what kind of ice strength we are using) … … 640 547 !! ** External : 641 548 !! 642 !! ** References : There are a lot of references and they are difficult to read 643 !! so the best thing to do is to accept what is here 644 !! 645 !! ** History : 646 !! This routine was freely inpired from CICE 647 !! authors William H. Lipscomb, LANL 648 !! Elizabeth C. Hunke (LANL) 649 !! the routines were available online 650 !! so, big thanks to the authors and to LANL 651 !! The modified features are the ridging/rafting algorithm 652 !! 653 !! (02-2006) Martin Vancoppenolle, UCL-ASTR 654 !! Voila, c'est la derniere etape du codage!!! J'en ai 655 !! plein le cul mais je m'en sors mieux qu'avant!!! 656 !!------------------------------------------------------------------ 549 !! ** References : 550 !! 551 !!---------------------------------------------------------------------- 657 552 !! * Arguments 658 553 … … 674 569 REAL(wp), DIMENSION(jpi,jpj) :: & 675 570 zworka !: temporary array used here 676 677 WRITE(numout,*) ' lim_itd_me_icestrength : compute ice strength '678 WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~ '679 571 680 572 !------------------------------------------------------------------------------! … … 707 599 hi * hi 708 600 709 ! strength(ji,jj) = strength(ji,jj) &710 ! + athorn(ji,jj,jl)/krdg(ji,jj,jl) &711 ! * 1.0/3.0 * (hrmax(ji,jj,jl)**3 - hrmin(ji,jj,jl)**3) &712 ! / (hrmax(ji,jj,jl)-hrmin(ji,jj,jl))713 601 !-------------------------- 714 602 ! PE gain from rafting ice … … 720 608 ! PE gain from ridging ice 721 609 !---------------------------- 722 ! anything from ridge porosity ?723 610 strength(ji,jj) = strength(ji,jj) & 724 611 + aridge(ji,jj,jl)/krdg(ji,jj,jl) & … … 753 640 END DO ! i 754 641 755 WRITE(numout,*) ' Pstar : ', Pstar756 WRITE(numout,*) ' vt_i : ', vt_i(jiindex,jjindex)757 WRITE(numout,*) ' C_rhg : ', C_rhg758 WRITE(numout,*) ' at_i : ', at_i(jiindex,jjindex)759 760 642 ksmooth = 1 761 643 … … 766 648 ! 5) Impact of brine volume 767 649 !------------------------------------------------------------------------------! 650 ! CAN BE REMOVED 768 651 ! 769 652 IF ( brinstren_swi .EQ. 1 ) THEN 770 771 WRITE(numout,*) ' mean brine volume : ', bv_i(jiindex,jjindex)772 653 773 654 DO jj = 1, jpj … … 778 659 zdummy = 0.0 779 660 ENDIF 780 ! zdummy = 0.0781 ! strength(ji,jj) = strength(ji,jj) * ( 1.0 - zdummy )782 ! Timco and O'brien, CRST, 1994, Flexural strength equation for sea ice783 661 strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bv_i(ji,jj),0.0))) 784 662 END DO ! j 785 663 END DO ! i 786 664 787 WRITE(numout,*) ' strength after brine volume correction : ', &788 strength(jiindex,jjindex)789 790 665 ENDIF 791 666 … … 861 736 ENDIF ! ksmooth 862 737 863 WRITE(numout,*) ' Ice strength : ', strength(jiindex,jjindex)864 865 738 ! Boundary conditions 866 739 CALL lbc_lnk( strength, 'T', 1. ) … … 885 758 !! ** External : 886 759 !! 887 !! ** History :888 !! This routine was freely inspired from CICE889 !! authors William H. Lipscomb, LANL890 !! Elizabeth C. Hunke (LANL)891 !! the routines were available online892 !! so, big thanks to the authors and to LANL893 !! The modified features are the ridging/rafting algorithm894 !!895 !! (02-2006) Martin Vancoppenolle, UCL-ASTR896 !! J'espere que nos lignes sont courbes897 !! Mais j'ai bien peur qu'elles ne soient droites898 760 !!---------------------------------------------------------------------! 899 761 !! * Arguments … … 901 763 INTEGER :: & 902 764 ji,jj, & ! horizontal indices 903 jl,jl1, & ! thickness category index 904 jm, & ! ice type index 905 index, & 765 jl, & ! thickness category index 906 766 krdg_index ! which participation function using 907 767 … … 913 773 Gsum ! Gsum(n) = sum of areas in categories 0 to n 914 774 915 REAL(wp), DIMENSION(jpi,jpj,-1:jpl) :: &916 Gsumra ! Gsum(n) = sum of areas in categories 0 to n917 918 REAL(wp), DIMENSION(jpi,jpj,jpm) :: &919 ai_typ ! Surface of respective ice types920 921 775 REAL(wp) :: & 922 776 hi, & ! ice thickness for each cat (m) … … 927 781 928 782 REAL(wp) :: & 929 zindb, zdummy, & 930 zweight, & 783 zdummy, & 931 784 epsi06 = 1.0e-6 932 785 933 INTEGER, DIMENSION(jpi,jpj,-1:jpl) :: &934 index_cumul ! temporary array used here935 786 !------------------------------------------------------------------------------! 936 787 … … 958 809 ! because of divergence during transport 959 810 960 !+++++961 ! 0.0 Compute surfaces of relative ice types962 ! BACK TO LIMA_MEC963 ! ai_typ(:,:,:) = 0.0964 ! DO jm = 1, jpm965 ! DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)966 ! DO jj = 1, jpj967 ! DO ji = 1, jpi968 ! ai_typ(ji,jj,jm) = ai_typ(ji,jj,jm) + a_i(ji,jj,jl)969 ! END DO970 ! END DO971 ! END DO972 ! END DO973 ! WRITE(numout,*) ' ai_typ ', ai_typ(jiindex,jjindex,1:3)974 ! BACK TO LIMA_MEC975 976 811 ! Compute cumulative thickness distribution function 977 812 ! Compute the cumulative thickness distribution function Gsum, 978 813 ! where Gsum(n) is the fractional area in categories 0 to n. 979 980 814 ! initial value (in h = 0) equals open water area 981 815 … … 1006 840 1007 841 ! Normalize the cumulative distribution to 1 1008 !former line1009 ! zworka(:,:) = 1.0 / Gsum(:,:,jpl)1010 WRITE(numout,*) ' AVANT : numit '1011 !new lines1012 842 DO jj = 1, jpj 1013 843 DO ji = 1, jpi 1014 ! IF (numit.GE.46365) THEN1015 IF ( (ji .EQ. jiindex) .AND. ( jj.EQ.jjindex ) ) THEN1016 WRITE(numout,*) ' ji ,jj : ',ji,jj, 'Gsum:', &1017 Gsum(ji,jj,jpl)1018 ENDIF1019 844 zworka(ji,jj) = 1.0 / Gsum(ji,jj,jpl) 1020 845 END DO 1021 846 END DO 1022 847 1023 WRITE(numout,*) ' APRES '1024 848 DO jl = 0, jpl 1025 849 DO jj = 1, jpj 1026 850 DO ji = 1, jpi 1027 IF ( (ji .EQ. jiindex) .AND. ( jj.EQ.jjindex ) ) THEN1028 WRITE(numout,*) ' ji ,jj : ',ji,jj,jl, 'Gsum:', &1029 Gsum(ji,jj,jl)1030 WRITE(numout,*) ' zworka ', zworka(ji,jj)1031 ENDIF1032 1033 851 Gsum(ji,jj,jl) = Gsum(ji,jj,jl) * zworka(ji,jj) 1034 852 END DO 1035 853 END DO 1036 854 END DO 1037 1038 ! BACK TO LIMA_MEC1039 ! For rafted ice1040 ! jm = 31041 ! DO jj = 1, jpj1042 ! DO ji = 1, jpi1043 ! Gsumra(ji,jj,ice_cat_bounds(jm,1)-1) = 0.01044 !! zworka(ji,jj) = 1.0 / ai_typ(ji,jj,jm) ! used to normalize Gsum1045 ! END DO1046 ! END DO1047 1048 ! DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)1049 ! DO jj = 1, jpj1050 ! DO ji = 1, jpi1051 ! Gsumra(ji,jj,jl) = Gsumra(ji,jj,jl-1) + a_i(ji,jj,jl)*zworka(ji,jj)1052 ! END DO1053 ! END DO1054 ! END DO1055 1056 ! +++++ [1057 WRITE(numout,*)1058 WRITE(numout,*) ' Gsum : ', Gsum(jiindex,jjindex,-1:jpl)1059 ! +++++ ]1060 1061 855 1062 856 ! 1.3 Compute participation function a(h) = b(h).g(h) (athorn) … … 1098 892 ! precompute exponential terms using Gsum as a work array 1099 893 zdummy = 1.0 / (1.0-EXP(-astari)) 1100 WRITE(numout,*) ' zdummy : ', zdummy1101 894 1102 895 DO jl = -1, jpl … … 1119 912 ENDIF ! krdg_index 1120 913 1121 WRITE(numout,*) ' astari : ', astari1122 WRITE(numout,*) ' athorn : ', athorn(jiindex,jjindex,0:jpl)1123 1124 ! BACK TO LIMA_MEC1125 ! jm = 31126 ! DO jl = ice_cat_bounds(3,1), ice_cat_bounds(3,2) ! only undeformed ice participates1127 ! DO jj = 1, jpj1128 ! DO ji = 1, jpi1129 ! IF (Gsumra(ji,jj,jl) < Gstar) THEN1130 ! athorn(ji,jj,jl) = Gstari * (Gsumra(ji,jj,jl)-Gsumra(ji,jj,jl-1)) * &1131 ! (2.0 - (Gsumra(ji,jj,jl-1)+Gsumra(ji,jj,jl))*Gstari)1132 ! ELSEIF (Gsumra(ji,jj,jl-1) < Gstar) THEN1133 ! athorn(ji,jj,jl) = Gstari * (Gstar-Gsumra(ji,jj,jl-1)) * &1134 ! (2.0 - (Gsumra(ji,jj,jl-1)+Gstar)*Gstari)1135 ! ELSE1136 ! athorn(ji,jj,jl) = 0.01137 ! ENDIF1138 ! END DO ! ji1139 ! END DO ! jj1140 ! END DO ! jl1141 ! BACK TO LIMA_MEC1142 1143 !!!1144 ! BACK TO LIMA_MEC1145 ! Weighting1146 ! Undeformed ice1147 ! DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2)1148 ! DO jj = 1, jpj1149 ! DO ji = 1, jpi1150 ! zindb = MAX( 0.0, SIGN( 1.0, a_i(ji,jj,jl) - epsi06 ) )1151 ! zweight = zindb*ai_typ(ji,jj,1) / &1152 ! MAX( ai_typ(ji,jj,1) + ai_typ(ji,jj,3), epsi06 )1153 ! athorn(ji,jj,jl) = zweight*athorn(ji,jj,jl)1154 ! END DO1155 ! END DO1156 ! END DO1157 1158 ! Rafted ice1159 ! jm = 31160 ! DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)1161 ! DO jj = 1, jpj1162 ! DO ji = 1, jpi1163 ! zindb = MAX( 0.0, SIGN( 1.0, a_i(ji,jj,jl) - epsi06 ) )1164 ! zweight = zindb*ai_typ(ji,jj,jm) / &1165 ! MAX( ai_typ(ji,jj,1) + ai_typ(ji,jj,3) , epsi06)1166 ! athorn(ji,jj,jl) = zweight*athorn(ji,jj,jl)1167 ! END DO1168 ! END DO1169 ! END DO1170 1171 !+++++1172 ! BACK TO LIMA_MEC1173 !!!1174 1175 914 ! Ridging and rafting ice participation functions 1176 915 IF ( raftswi .EQ. 1 ) THEN … … 1224 963 1225 964 ENDIF 1226 !+++++1227 WRITE(numout,*)1228 WRITE(numout,*) ' aridge : ', aridge(jiindex, jjindex, 1:jpl)1229 WRITE(numout,*) ' araft : ', araft (jiindex, jjindex, 1:jpl)1230 !+++++1231 965 1232 966 !----------------------------------------------------------------- … … 1290 1024 DO ji = 1, jpi 1291 1025 aksum(ji,jj) = aksum(ji,jj) & 1292 ! + athorn(ji,jj,jl) * (1.0 - 1.0/krdg(ji,jj,jl))1293 1026 + aridge(ji,jj,jl) * (1.0 - 1.0/krdg(ji,jj,jl)) & 1294 1027 + araft (ji,jj,jl) * (1.0 - 1.0/kraft) … … 1297 1030 END DO 1298 1031 1299 WRITE(numout,*) ' krdg : ', krdg(jiindex,jjindex,1:jpl)1300 WRITE(numout,*) ' hraft : ', hraft(jiindex,jjindex,1:jpl)1301 WRITE(numout,*) ' aksum : ', aksum(jiindex,jjindex)1302 1303 1032 END SUBROUTINE lim_itd_me_ridgeprep 1304 1033 … … 1324 1053 !! ** External : 1325 1054 !! 1326 !! ** History :1327 !! This routine was freely inpired from CICE1328 !! authors William H. Lipscomb, LANL1329 !! Elizabeth C. Hunke (LANL)1330 !! the routines were available online1331 !! so, big thanks to the authors and to LANL1332 !! The modified features are the ridging/rafting algorithm1333 !!1334 !! (02-2006) Martin Vancoppenolle, UCL-ASTR1335 !! Dans 5 mois je suis sensé avoir fini ma thèse alors j'ai interet a cavaler1336 !! sinon ca va mal se passer1337 !! J'ai trop chaud, je suis fatigué1338 !! Mon ame se desseche t elle ??1339 !! Suis je en train de changer ???1340 !! Oublie-je l'essentiel de mon histoire ??1341 !!1342 1055 1343 1056 REAL (wp), DIMENSION(jpi,jpj), INTENT(IN) :: & … … 1355 1068 jk, & ! ice layer index 1356 1069 ij, & ! horizontal index, combines i and j loops 1357 icells, & ! number of cells with aicen > puny 1358 index ! index for debug ice points 1070 icells ! number of cells with aicen > puny 1359 1071 1360 1072 INTEGER, DIMENSION(1:(jpi+1)*(jpj+1)) :: & … … 1370 1082 vsnon_init, & ! snow volume before ridging 1371 1083 esnon_init, & ! snow energy before ridging 1372 smv_i_init, 1373 oa_i_init 1084 smv_i_init, & ! ice salinity before ridging 1085 oa_i_init ! ice age before ridging 1374 1086 1375 1087 REAL(wp), DIMENSION(jpi,jpj,jkmax,jpl) :: & … … 1380 1092 ardg1 , & ! area of ice ridged 1381 1093 ardg2 , & ! area of new ridges 1382 virdg , & ! ice volume of ridging ice1383 1094 vsrdg , & ! snow volume of ridging ice 1384 1095 esrdg , & ! snow energy of ridging ice 1385 smrdg , & ! salinity of ridging ice1386 oirdg , & ! volumic age content of ridging ice1387 1096 oirdg1 , & ! areal age content of ridged ice 1388 1097 oirdg2 , & ! areal age content of ridging ice … … 1397 1106 srdg1 , & ! sal*volume of ice ridged 1398 1107 srdg2 , & ! sal*volume of new ridges 1399 smsw , & ! sal*volume of water trapped into ridges 1400 totvrdg1 , & ! total volume of ice ridged 1401 totvrdg2 ! total volume of new ridges 1108 smsw ! sal*volume of water trapped into ridges 1402 1109 1403 1110 REAL(wp), DIMENSION(jpi,jpj) :: & … … 1409 1116 esrft , & ! snow energy of rafting ice 1410 1117 smrft , & ! salinity of rafting ice 1411 oirft , & ! volumic age content of rafting ice1412 1118 oirft1 , & ! areal age content of rafted ice 1413 1119 oirft2 ! areal age content of rafting ice 1414 1120 1415 1121 REAL(wp), DIMENSION(jpi,jpj,jkmax) :: & 1416 eirdg , & ! ice energy of ridging ice1417 1122 eirft , & ! ice energy of rafting ice 1418 1123 erdg1 , & ! enth*volume of ice ridged … … 1477 1182 END DO !jj 1478 1183 END DO !ji 1479 1480 !+++++1481 ! WRITE(numout,*) ' 1. Change in open water area '1482 ! WRITE(numout,*) ' closing_gross', closing_gross(jiindex, jjindex)1483 ! WRITE(numout,*) ' opning ', opning (jiindex, jjindex)1484 ! WRITE(numout,*) ' ato_i ', ato_i (jiindex, jjindex)1485 ! WRITE(numout,*)1486 !+++++1487 1184 1488 1185 ! if negative open water area alert it … … 1527 1224 END DO !jl 1528 1225 1529 !+++++ 1530 WRITE(numout,*) ' State variables before ridging is computed ' 1531 WRITE(numout,*) ' a_i ', a_i(jiindex, jjindex, 1:jpl) 1532 WRITE(numout,*) ' v_i ', v_i(jiindex, jjindex, 1:jpl) 1533 WRITE(numout,*) ' v_s ', v_s(jiindex, jjindex, 1:jpl) 1534 WRITE(numout,*) ' e_s ', e_s(jiindex, jjindex, 1, 1:jpl) 1535 WRITE(numout,*) ' smv_i ', smv_i(jiindex, jjindex, 1:jpl) 1536 WRITE(numout,*) ' oa_i ', oa_i (jiindex, jjindex, 1:jpl) 1537 WRITE(numout,*) ' e_i 1 ', e_i (jiindex, jjindex, 1, 1:jpl) 1538 WRITE(numout,*) ' e_i 2 ', e_i (jiindex, jjindex, 2, 1:jpl) 1539 WRITE(numout,*) 1540 ! !+++++ 1541 1226 ! 1542 1227 !----------------------------------------------------------------- 1543 1228 ! 3) Pump everything from ice which is being ridged / rafted … … 1552 1237 !------------------------------------------------ 1553 1238 1554 !+++++1555 ! WRITE(numout,*)1556 ! WRITE(numout,*) ' *** jl1 *** ', jl11557 !+++++1558 1239 icells = 0 1559 1240 DO jj = 1, jpj … … 1571 1252 large_afrft = .false. 1572 1253 1573 ! !+++++1574 ! WRITE(numout,*) ' 3.1 ridging cells '1575 ! WRITE(numout,*) ' icells ', icells1576 ! WRITE(numout,*)1577 !+++++1578 1579 1254 DO ij = 1, icells 1580 1255 ji = indxi(ij) … … 1594 1269 oirdg2(ji,jj)= oirdg1(ji,jj) / krdg(ji,jj,jl1) 1595 1270 oirft2(ji,jj)= oirft1(ji,jj) / kraft 1596 1597 !+++++1598 ! IF ( ( ji.EQ.jiindex ) .AND. (jj.EQ.jjindex) ) THEN1599 ! WRITE(numout,*) ' 3.2 area of ridging ice and of new ridge '1600 ! WRITE(numout,*) ' aridge ', aridge(ji,jj,jl1)1601 ! WRITE(numout,*) ' araft ', araft (ji,jj,jl1)1602 ! WRITE(numout,*) ' krdg ', krdg (ji,jj,jl1)1603 ! WRITE(numout,*) ' kraft ', kraft1604 ! WRITE(numout,*) ' clogross ', closing_gross(ji,jj)1605 ! WRITE(numout,*) ' ardg1 : ', ardg1(jiindex,jjindex)1606 ! WRITE(numout,*) ' arft1 : ', arft1(jiindex,jjindex)1607 ! WRITE(numout,*) ' ardg2 : ', ardg2(jiindex,jjindex)1608 ! WRITE(numout,*) ' arft2 : ', arft2(jiindex,jjindex)1609 ! WRITE(numout,*)1610 ! ENDIF1611 !+++++1612 1271 1613 1272 !--------------------------------------------------------------- … … 1629 1288 ENDIF 1630 1289 1631 ! !+++++1632 ! IF ( ( ji.EQ.jiindex ) .AND. (jj.EQ.jjindex) ) THEN1633 ! WRITE(numout,*) ' 3.3 Ridging / rafting fractions '1634 ! WRITE(numout,*) ' afrac : ', afrac(ji,jj)1635 ! WRITE(numout,*) ' afrft : ', afrft(ji,jj)1636 ! WRITE(numout,*)1637 ! ENDIF1638 ! !+++++1639 1640 1290 !-------------------------------------------------------------------------- 1641 1291 ! 3.4) Subtract area, volume, and energy from ridging 1642 1292 ! / rafting category n1. 1643 1293 !-------------------------------------------------------------------------- 1644 ! jl1 looping 1-jpl1645 ! ij looping 1-icells1646 1647 ! ridging volumes, heat contents ...1648 ! virdg is the volume of ridged ice, at the end of ridging process1649 ! i.e., it includes 30% porosity1650 ! virdg(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj)1651 1294 vrdg1(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj) / & 1652 1295 ( 1.0 + ridge_por ) … … 1659 1302 ( 1. + ridge_por ) 1660 1303 srdg2(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) 1661 ! oirdg(ji,jj) = oa_i_init(ji,jj,jl1) * afrac(ji,jj) * &1662 ! ( 1.0 - ridge_por )1663 1304 1664 1305 ! rafting volumes, heat contents ... … … 1667 1308 esrft(ji,jj) = esnon_init(ji,jj,jl1) * afrft(ji,jj) 1668 1309 smrft(ji,jj) = smv_i_init(ji,jj,jl1) * afrft(ji,jj) 1669 ! oirft(ji,jj) = oa_i_init(ji,jj,jl1) * afrft(ji,jj)1670 1310 1671 1311 ! substract everything … … 1683 1323 ! Salinity 1684 1324 !------------- 1685 ! salinity times volume of the seawater into ridges1686 ! WRITE(numout,*) ' sss_io : ', sss_io(ji,jj)1687 ! WRITE(numout,*) ' vsw : ', vsw(ji,jj)1688 ! WRITE(numout,*) ' ridge_por : ', ridge_por1689 1690 1325 smsw(ji,jj) = sss_io(ji,jj) * vsw(ji,jj) * ridge_por 1691 1326 … … 1704 1339 srdg2(ji,jj) = sm_newridge * vrdg2(ji,jj) 1705 1340 1706 1707 1341 !------------------------------------ 1708 ! 3. 5Increment ridging diagnostics1342 ! 3.6 Increment ridging diagnostics 1709 1343 !------------------------------------ 1710 1344 … … 1720 1354 1721 1355 !------------------------------------------ 1722 ! 3. 6Put the snow somewhere in the ocean1356 ! 3.7 Put the snow somewhere in the ocean 1723 1357 !------------------------------------------ 1724 1358 … … 1742 1376 1743 1377 !----------------------------------------------------------------- 1744 ! 3. 7)Compute quantities used to apportion ice among categories1378 ! 3.8 Compute quantities used to apportion ice among categories 1745 1379 ! in the n2 loop below 1746 1380 !----------------------------------------------------------------- … … 1754 1388 1755 1389 1756 !++++++++1757 IF ( vrdg1(ji,jj) .GT. 0.0 ) THEN1758 IF ( ( vrdg2(ji,jj) / vrdg1(ji,jj) - ( 1. + ridge_por ) .GT. epsi10 ) ) &1759 THEN1760 WRITE(numout,*) ' ALERTE in ridging routine '1761 WRITE(numout,*) ' difference : ', vrdg2(ji,jj) / vrdg1(ji,jj) - ( 1. + ridge_por )1762 WRITE(numout,*) ' vrdg2/vrdg1 ', vrdg2(ji,jj) / &1763 vrdg1(ji,jj)1764 WRITE(numout,*) ' should be equal to 1+p : ', 1 + ridge_por1765 WRITE(numout,*) ' vrdg1 : ', vrdg1(ji,jj)1766 WRITE(numout,*) ' vrdg2 : ', vrdg2(ji,jj)1767 ENDIF1768 ENDIF1769 1770 !++++++++1771 !1772 1773 1390 END DO ! ij 1774 1391 1775 ! !+++++1776 ! !+++++1777 1778 ! !+++++1779 ! WRITE(numout,*)1780 ! WRITE(numout,*) ' Mass of salt of ridging ice 1 : ', srdg1(jiindex,jjindex)1781 ! WRITE(numout,*) ' Mass of salt of ridging ice 2 : ', srdg2(jiindex,jjindex)1782 ! WRITE(numout,*) ' Salinity of ridging ice : ', smrdg(jiindex,jjindex) / &1783 ! virdg(jiindex,jjindex) / (1.0-ridge_por)1784 ! WRITE(numout,*) ' Mass of salt after ridging : ', sm_i(ji,jj,jl1)1785 ! WRITE(numout,*) ' Sal after removal of ridg.ice ', sm_i(ji,jj,jl1) / &1786 ! vicen_init(ji,jj,jl1)1787 !+++++1788 1789 1790 1392 !-------------------------------------------------------------------- 1791 ! 3. 8)Compute ridging ice enthalpy, remove it from ridging ice and1393 ! 3.9 Compute ridging ice enthalpy, remove it from ridging ice and 1792 1394 ! compute ridged ice enthalpy 1793 1395 !-------------------------------------------------------------------- … … 1838 1440 ji = indxi(ij) 1839 1441 jj = indxj(ij) 1840 ! ridging diagnostics1841 ! WRITE(numout,*) ' debut ', eice_init(ji,jj)1842 ! WRITE(numout,*) ' erdg2 ', erdg2(ji,jj,jk)1843 ! WRITE(numout,*) ' erdg1 ', erdg1(ji,jj,jk)1844 ! WRITE(numout,*) ' ersw ', ersw(ji,jj,jk)1845 ! WRITE(numout,*) ' vsw ', vsw(ji,jj)1846 1442 eice_init(ji,jj) = eice_init(ji,jj) + erdg2(ji,jj,jk) - & 1847 1443 erdg1(ji,jj,jk) … … 1849 1445 END DO !jk 1850 1446 ENDIF 1851 1852 !+++++ [1853 ! WRITE(numout,*)1854 ! WRITE(numout,*) ' ridging ice enthalpy ', jl11855 ! WRITE(numout,*) ' eird1 ', erdg1(jiindex,jjindex,1:nlay_i)1856 ! WRITE(numout,*) ' eird2 old ', erdg2(jiindex,jjindex,1:nlay_i)1857 ! WRITE(numout,*) ' eirft ', eirft(jiindex,jjindex,1:2)1858 ! WRITE(numout,*) ' e_i ', e_i (jiindex,jjindex,1:nlay_i,jl1)1859 !+++++ ]1860 1447 1861 1448 IF (large_afrac) THEN ! there is a bug … … 1884 1471 ENDIF ! large_afrft 1885 1472 1886 !+++++ [1887 ! IF ( ( ji.EQ.jiindex ) .AND. (jj.EQ.jjindex) ) THEN1888 ! WRITE(numout,*) ' 3.4 Substract area, volume and energy '1889 ! WRITE(numout,*) ' virdg : ', virdg(ji,jj)1890 ! WRITE(numout,*) ' vsrdg ', vsrdg(ji,jj)1891 ! WRITE(numout,*) ' esrdg ', esrdg(ji,jj)1892 ! WRITE(numout,*) ' smrdg ', smrdg(ji,jj)1893 ! WRITE(numout,*) ' oirdg ', oirdg(ji,jj)1894 ! WRITE(numout,*)1895 ! WRITE(numout,*) ' virft : ', virft(ji,jj)1896 ! WRITE(numout,*) ' vsrft ', vsrft(ji,jj)1897 ! WRITE(numout,*) ' esrft ', esrft(ji,jj)1898 ! WRITE(numout,*) ' smrft ', smrft(ji,jj)1899 ! WRITE(numout,*) ' oirft ', oirft(ji,jj)1900 ! WRITE(numout,*)1901 1902 ! WRITE(numout,*) ' a_i : ', a_i(ji,jj,jl1)1903 ! WRITE(numout,*) ' v_i : ', v_i(ji,jj,jl1)1904 ! WRITE(numout,*) ' v_s ', v_s(ji,jj,jl1)1905 ! WRITE(numout,*) ' e_s ', e_s(ji,jj,1,jl1)1906 ! WRITE(numout,*) ' smv_i ', smv_i(ji,jj,jl1)1907 ! WRITE(numout,*) ' oa_i ', oa_i(ji,jj,jl1)1908 ! ENDIF1909 !+++++ ]1910 1911 1473 !------------------------------------------------------------------------------- 1912 1474 ! 4) Add area, volume, and energy of new ridge to each category jl2 … … 1946 1508 1947 1509 END DO ! ij 1948 ! WRITE(numout,*) ' new ridges being built from ', jl11949 ! WRITE(numout,*) ' added to ', jl21950 ! WRITE(numout,*) ' a_i : ', a_i(jiindex,jjindex,jl2)1951 ! WRITE(numout,*) ' v_i : ', v_i(jiindex,jjindex,jl2)1952 1510 1953 1511 ! Transfer ice energy to category jl2 by ridging … … 1964 1522 END DO ! jl2 (new ridges) 1965 1523 1966 ! WRITE(numout,*) ' new ridges built '1967 ! WRITE(numout,*) ' a_i : ', a_i(jiindex,jjindex,jl1)1968 ! WRITE(numout,*) ' v_i : ', v_i(jiindex,jjindex,jl1)1969 ! WRITE(numout,*) ' e_i : ', e_i(jiindex,jjindex,:,jl1)1970 1971 1524 DO jl2 = ice_cat_bounds(1,1), ice_cat_bounds(1,2) 1972 ! jl2 = ice_cat_bounds(3,1) ! rafted category1973 1525 1974 1526 DO ij = 1, icells … … 1990 1542 oa_i(ji,jj,jl2) = oa_i(ji,jj,jl2) & 1991 1543 + oirft2(ji,jj) 1992 1993 ! IF ( (ji.EQ.jiindex) .AND. (jj.EQ.jjindex) ) THEN1994 ! WRITE(numout,*) ' raft ice from jl1 ', jl1 , &1995 ! ' to jl2 :', jl21996 ! WRITE(numout,*) ' hraft (jl1) ', hraft(ji,jj,jl1)1997 ! WRITE(numout,*) ' hi_maxs ', hi_max(jl2-1), hi_max(jl2)1998 ! ENDIF1999 2000 1544 ENDIF ! hraft 2001 1545 … … 2017 1561 END DO ! jl2 2018 1562 2019 ! WRITE(numout,*) ' new rafted ice built '2020 ! WRITE(numout,*) ' a_i : ', a_i(jiindex,jjindex,jl1)2021 ! WRITE(numout,*) ' v_i : ', v_i(jiindex,jjindex,jl1)2022 ! WRITE(numout,*) ' e_i : ', e_i(jiindex,jjindex,:,jl1)2023 2024 1563 END DO ! jl1 (deforming categories) 2025 1564 … … 2060 1599 2061 1600 INTEGER :: ji, jj, jl 2062 2063 ! WRITE(numout,*) ' lim_itd_me_asumr : sum the ice / open water fractions'2064 ! WRITE(numout,*) ' ~~~~~~~~~~~~~~~~ '2065 1601 2066 1602 !----------------------------------------------------------------- … … 2172 1708 xtmp ! temporary variable 2173 1709 2174 WRITE(numout,*) ' lim_itd_me_zapsmall '2175 WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~ '2176 2177 1710 DO jl = 1, jpl 2178 1711
Note: See TracChangeset
for help on using the changeset viewer.
