- Timestamp:
- 2016-07-01T18:02:45+02:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
r5602 r6772 28 28 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 29 29 USE wrk_nemo ! work arrays 30 USE fldread ! read input fields 31 USE iom 30 32 31 33 IMPLICIT NONE … … 47 49 REAL(wp) :: rn_tmi_ini_s ! initial temperature 48 50 51 INTEGER , PARAMETER :: jpfldi = 7 ! maximum number of files to read 52 INTEGER , PARAMETER :: jp_hicif = 1 ! index of thick (m) at T-point 53 INTEGER , PARAMETER :: jp_hsnif = 2 ! index of thick (m) at T-point 54 INTEGER , PARAMETER :: jp_frld = 3 ! index of ice fraction (%) at T-point 55 INTEGER , PARAMETER :: jp_sist = 4 ! index of ice surface temp (K) at T-point 56 INTEGER , PARAMETER :: jp_tbif1 = 5 ! index of ice temp lev1 (K) at T-point 57 INTEGER , PARAMETER :: jp_tbif2 = 6 ! index of ice temp lev2 (K) at T-point 58 INTEGER , PARAMETER :: jp_tbif3 = 7 ! index of ice temp lev3 (K) at T-point 59 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: si ! structure of input fields (file informations, fields read) 60 61 REAL(wp),DIMENSION(:,:) ,ALLOCATABLE :: hicif_ini,hsnif_ini,frld_ini,sist_ini, zswitch 62 REAL(wp),DIMENSION(:,:,:),ALLOCATABLE :: tbif_ini 63 49 64 LOGICAL :: ln_iceini ! initialization or not 65 LOGICAL :: ln_limini_file ! Ice initialization state from 2D netcdf file 50 66 !!---------------------------------------------------------------------- 51 67 !! LIM 3.0, UCL-LOCEAN-IPSL (2008) … … 91 107 REAL(wp), POINTER, DIMENSION(:) :: zht_i_ini, zat_i_ini, zvt_i_ini, zht_s_ini, zsm_i_ini, ztm_i_ini 92 108 REAL(wp), POINTER, DIMENSION(:,:) :: zh_i_ini, za_i_ini, zv_i_ini 93 REAL(wp), POINTER, DIMENSION(:,:) :: zswitch ! ice indicator94 109 INTEGER, POINTER, DIMENSION(:,:) :: zhemis ! hemispheric index 95 110 !-------------------------------------------------------------------- 96 111 97 CALL wrk_alloc( jpi, jpj, zswitch )98 112 CALL wrk_alloc( jpi, jpj, zhemis ) 99 113 CALL wrk_alloc( jpl, 2, zh_i_ini, za_i_ini, zv_i_ini ) … … 150 164 ! 3) Initialization of sea ice state variables 151 165 !-------------------------------------------------------------------- 166 IF( ln_limini_file )THEN 167 168 CALL limini_file 169 170 ELSE 152 171 153 172 !----------------------------- … … 376 395 tn_ice (:,:,:) = t_su (:,:,:) 377 396 397 ENDIF !ln_limini_file 398 378 399 ELSE 379 400 ! if ln_iceini=false … … 399 420 END DO 400 421 END DO 401 422 402 423 ENDIF ! ln_iceini 403 424 … … 451 472 452 473 453 CALL wrk_dealloc( jpi, jpj, zswitch )454 474 CALL wrk_dealloc( jpi, jpj, zhemis ) 455 475 CALL wrk_dealloc( jpl, 2, zh_i_ini, za_i_ini, zv_i_ini ) … … 474 494 !! 8.5 ! 07-11 (M. Vancoppenolle) rewritten initialization 475 495 !!----------------------------------------------------------------------------- 476 NAMELIST/namiceini/ ln_iceini, rn_thres_sst, rn_hts_ini_n, rn_hts_ini_s, rn_hti_ini_n, rn_hti_ini_s, & 477 & rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, rn_smi_ini_s, rn_tmi_ini_n, rn_tmi_ini_s 478 INTEGER :: ios ! Local integer output status for namelist read 496 ! 497 INTEGER :: ios,ierr,inum_ice ! Local integer output status for namelist read 498 INTEGER :: ji,jj 499 INTEGER :: ifpr, ierror 500 ! 501 CHARACTER(len=100) :: cn_dir ! Root directory for location of ice files 502 TYPE(FLD_N) :: sn_hicif, sn_hsnif, sn_frld, sn_sist 503 TYPE(FLD_N) :: sn_tbif1, sn_tbif2, sn_tbif3 504 TYPE(FLD_N), DIMENSION(jpfldi) :: slf_i ! array of namelist informations on the fields to read 505 ! 506 NAMELIST/namiceini/ ln_iceini, ln_limini_file, rn_thres_sst, rn_hts_ini_n, rn_hts_ini_s, & 507 & rn_hti_ini_n, rn_hti_ini_s, rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, & 508 & rn_smi_ini_s, rn_tmi_ini_n, rn_tmi_ini_s, & 509 & sn_hicif, sn_hsnif, sn_frld, sn_sist, & 510 & sn_tbif1, sn_tbif2, sn_tbif3, cn_dir 479 511 !!----------------------------------------------------------------------------- 480 512 ! … … 488 520 IF(lwm) WRITE ( numoni, namiceini ) 489 521 522 slf_i(jp_hicif) = sn_hicif ; slf_i(jp_hsnif) = sn_hsnif 523 slf_i(jp_frld) = sn_frld ; slf_i(jp_sist) = sn_sist 524 slf_i(jp_tbif1) = sn_tbif1 ; slf_i(jp_tbif2) = sn_tbif2 ; slf_i(jp_tbif3) = sn_tbif3 525 490 526 ! Define the initial parameters 491 527 ! ------------------------- … … 496 532 WRITE(numout,*) '~~~~~~~~~~~~~~~' 497 533 WRITE(numout,*) ' initialization with ice (T) or not (F) ln_iceini = ', ln_iceini 534 WRITE(numout,*) ' initialization with ice (T) or not (F) ln_limini_file = ', ln_limini_file 498 535 WRITE(numout,*) ' threshold water temp. for initial sea-ice rn_thres_sst = ', rn_thres_sst 499 536 WRITE(numout,*) ' initial snow thickness in the north rn_hts_ini_n = ', rn_hts_ini_n … … 509 546 ENDIF 510 547 548 IF( ln_limini_file ) THEN ! Ice initialization using input file 549 ! 550 ierr = alloc_lim_istate_init() 551 ! 552 ! CALL iom_open( 'Ice_initialization.nc', inum_ice ) 553 ! ! 554 ! IF( inum_ice > 0 ) THEN 555 ! IF(lwp) WRITE(numout,*) 556 ! IF(lwp) WRITE(numout,*) ' ice state initialization with : Ice_initialization.nc' 557 ! 558 ! CALL iom_get( inum_ice, jpdom_data, 'hicif', hicif_ini ) 559 ! CALL iom_get( inum_ice, jpdom_data, 'hsnif', hsnif_ini ) 560 ! CALL iom_get( inum_ice, jpdom_data, 'frld' , frld_ini ) 561 ! CALL iom_get( inum_ice, jpdom_data, 'ts' , sist_ini ) 562 ! CALL iom_get( inum_ice, jpdom_unknown, 'tbif', tbif_ini(1:nlci,1:nlcj,:), & 563 ! & kstart = (/ mig(1),mjg(1),1 /), kcount = (/ nlci,nlcj,3 /) ) 564 ! ! put some values in the extra-halo... 565 566 ! set si structure 567 ALLOCATE( si(jpfldi), STAT=ierror ) 568 IF( ierror > 0 ) THEN 569 CALL ctl_stop( 'Ice_ini in limistate: unable to allocate si structure' ) ; RETURN 570 ENDIF 571 572 DO ifpr= 1, jpfldi 573 ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) ) 574 ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) ) 575 END DO 576 577 ! fill si with slf_i and control print 578 CALL fld_fill( si, slf_i, cn_dir, 'lim_istate', 'lim istate ini', 'numnam_ice' ) 579 580 CALL fld_read( nit000, 1, si ) ! input fields provided at the current time-step 581 582 hicif_ini(:,:) = si(jp_hicif)%fnow(:,:,1) 583 hsnif_ini(:,:) = si(jp_hsnif)%fnow(:,:,1) 584 frld_ini(:,:) = si(jp_frld)%fnow(:,:,1) 585 sist_ini(:,:) = si(jp_sist)%fnow(:,:,1) 586 tbif_ini(:,:,1) = si(jp_tbif1)%fnow(:,:,1) 587 tbif_ini(:,:,2) = si(jp_tbif2)%fnow(:,:,1) 588 tbif_ini(:,:,3) = si(jp_tbif3)%fnow(:,:,1) 589 590 DO jj = nlcj+1, jpj ; tbif_ini(1:nlci,jj,:) = tbif_ini(1:nlci,nlej,:) ; END DO 591 DO ji = nlci+1, jpi ; tbif_ini(ji ,: ,:) = tbif_ini(nlei ,: ,:) ; END DO 592 593 ! CALL iom_close( inum_ice) 594 ! ! 595 ! ENDIF 596 ENDIF 597 511 598 END SUBROUTINE lim_istate_init 512 599 600 SUBROUTINE limini_file 601 !!----------------------------------------------------------------------------- 602 !! 603 !! 604 !! 605 !! 606 !!----------------------------------------------------------------------------- 607 INTEGER :: jl,ji,jj,jk 608 INTEGER :: jl0 609 INTEGER :: i_fill,jit,jjt 610 REAL(wp) :: ztest_1, ztest_2, ztest_3, ztest_4, ztests, zsigma, zarg, zA, zV, zA_cons, zV_cons, zconv,zH 611 REAL(wp) :: eps=1.e-6 612 REAL(wp) :: zmin,zmax 613 !rbb REAL(wp) :: epsi20,ztmelts,zdh 614 REAL(wp) ::ztmelts,zdh 615 616 REAL(wp), POINTER, DIMENSION(:,:) :: zhm_i_ini, zat_i_ini, zvt_i_ini, zhm_s_ini, zsm_i_ini 617 REAL(wp), POINTER, DIMENSION(:,:,:) :: zv_i_ini 618 REAL(wp), POINTER, DIMENSION(:,:,:) :: zht_i_ini,za_i_ini 619 REAL(wp), POINTER, DIMENSION(:,:) :: zidto ! ice indicator 620 !----------------------------------------------------------------------------- 621 IF(lwp)WRITE(numout,*)"limistate: read file : " 622 623 CALL wrk_alloc(jpl,jpi,jpj, zv_i_ini) 624 CALL wrk_alloc( jpi,jpj, zhm_i_ini, zat_i_ini, zvt_i_ini, zhm_s_ini, zsm_i_ini ) 625 CALL wrk_alloc( jpl,jpi,jpj,zht_i_ini,za_i_ini) 626 CALL wrk_alloc( jpi,jpj,zidto ) 627 628 zhm_i_ini(:,:) = hicif_ini(:,:) ! ice thickness 629 zat_i_ini(:,:) = 1._wp - frld_ini(:,:) ! ice concentration 630 zvt_i_ini(:,:) = zhm_i_ini(:,:) * zat_i_ini(:,:) ! ice volume 631 zhm_s_ini(:,:) = hsnif_ini(:,:) ! snow depth 632 633 zht_i_ini(:,:,:) = 0._wp 634 za_i_ini(:,:,:) = 0._wp 635 zv_i_ini(:,:,:) = 0._wp 636 637 zat_i_ini(:,:) = MIN( zat_i_ini(:,:) , 1.0_wp ) 638 639 640 DO ji = 1, jpi 641 DO jj = 1, jpj 642 643 IF( zat_i_ini(ji,jj) .GT. 0._wp .AND. zhm_i_ini(ji,jj) .GT. 0._wp )THEN 644 645 646 IF( gphit(ji,jj) .GE. 0._wp )THEN ; zsm_i_ini(ji,jj) = rn_smi_ini_n 647 ELSE ; zsm_i_ini(ji,jj) = rn_smi_ini_s 648 ENDIF 649 650 jl0 = 1 651 DO jl = 2, jpl 652 IF ( ( zhm_i_ini(ji,jj) .GT. hi_max(jl-1) ) .AND. & 653 ( zhm_i_ini(ji,jj) .LE. hi_max(jl) ) ) THEN 654 jl0 = jl 655 ENDIF 656 END DO 657 658 IF( jl0==1 )THEN 659 660 zht_i_ini(1,ji,jj) = zhm_i_ini(ji,jj) 661 za_i_ini(1,ji,jj) = zat_i_ini(ji,jj) 662 zht_i_ini(2:jpl,ji,jj) = 0._wp 663 za_i_ini(2:jpl,ji,jj) = 0._wp 664 665 ELSE ! jl0 ne 1 666 ztest_1 = 0 ; ztest_2 = 0 ; ztest_3 = 0 ; ztest_4 = 0 667 668 DO i_fill = jpl, 1, -1 669 IF( ( ztest_1 + ztest_2 + ztest_3 + ztest_4 ) .NE. 4 ) THEN 670 671 !---------------------------- 672 ! fill the i_fill categories 673 !---------------------------- 674 ! *** 1 category to fill 675 IF( i_fill .EQ. 1 ) THEN 676 zht_i_ini(1,ji,jj) = zhm_i_ini(ji,jj) 677 za_i_ini(1,ji,jj) = zat_i_ini(ji,jj) 678 zht_i_ini(2:jpl,ji,jj) = 0._wp 679 za_i_ini(2:jpl,ji,jj) = 0._wp 680 ELSE 681 682 ! *** >1 categores to fill 683 !--- Ice thicknesses in the i_fill - 1 first categories 684 DO jl = 1, i_fill - 1 685 zht_i_ini(jl,ji,jj) = 0.5 * ( hi_max(jl) + hi_max(jl-1) ) 686 END DO 687 688 !--- jl0: most likely index where cc will be maximum 689 DO jl = 1, jpl 690 IF ( ( zhm_i_ini(ji,jj) .GT. hi_max(jl-1) ) .AND. & 691 ( zhm_i_ini(ji,jj) .LE. hi_max(jl) ) ) THEN 692 jl0 = jl 693 ENDIF 694 END DO 695 jl0 = MIN(jl0, i_fill) 696 697 !--- Concentrations 698 za_i_ini(jl0,ji,jj) = zat_i_ini(ji,jj) / SQRT(REAL(jpl)) 699 DO jl = 1, i_fill - 1 700 IF ( jl .NE. jl0 ) THEN 701 zsigma = 0.5 * zhm_i_ini(ji,jj) 702 zarg = ( zht_i_ini(jl,ji,jj) - zhm_i_ini(ji,jj) ) / zsigma 703 za_i_ini(jl,ji,jj) = za_i_ini(jl0,ji,jj) * EXP(-zarg**2) 704 ENDIF 705 END DO 706 707 zA = 0. ! sum of the areas in the jpl categories 708 DO jl = 1, i_fill - 1 709 zA = zA + za_i_ini(jl,ji,jj) 710 END DO 711 za_i_ini(i_fill,ji,jj) = zat_i_ini(ji,jj) - zA ! ice conc in the last category 712 IF ( i_fill .LT. jpl ) za_i_ini(i_fill+1:jpl, ji,jj) = 0._wp 713 714 !--- Ice thickness in the last category 715 zV = 0. ! sum of the volumes of the N-1 categories 716 DO jl = 1, i_fill - 1 717 zV = zV + za_i_ini(jl,ji,jj)*zht_i_ini(jl,ji,jj) 718 END DO 719 zht_i_ini(i_fill,ji,jj) = ( zvt_i_ini(ji,jj) - zV ) /za_i_ini(i_fill,ji,jj) 720 IF ( i_fill .LT. jpl ) zht_i_ini(i_fill+1:jpl, ji,jj) = 0._wp 721 722 !--- volumes 723 zv_i_ini(:,ji,jj) = za_i_ini(:,ji,jj) * zht_i_ini(:,ji,jj) 724 IF ( i_fill .LT. jpl ) zv_i_ini(i_fill+1:jpl, ji,jj) = 0._wp 725 726 ENDIF ! i_fill 727 728 !--------------------- 729 ! Compatibility tests 730 !--------------------- 731 ! Test 1: area conservation 732 zA_cons = SUM(za_i_ini(:,ji,jj)) ; zconv = ABS(zat_i_ini(ji,jj) - zA_cons ) 733 IF ( zconv .LT. 1.0e-6 ) THEN 734 ztest_1 = 1 735 ELSE 736 ! this write is useful 737 !WRITE(numout,*) ' * TEST1 AREA NOT CONSERVED *** zA_cons = ', zA_cons,' zat_i_ini = ',zat_i_ini(ji,jj) 738 !WRITE(numout,*) 'ji,jj,narea ',ji,jj,narea 739 !WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj) 740 !WRITE(numout,*) ' zhm_i_ini : ', zhm_i_ini(ji,jj) 741 !WRITE(numout,*) ' zht_i_ini(:,jij,jj) ',zht_i_ini(:,ji,jj) 742 !WRITE(numout,*) ' za_i_ini(:,jij,jj) ',za_i_ini(:,ji,jj) 743 !WRITE(numout,*) ' hi_max ',hi_max 744 !WRITE(numout,*) ' jl0 = ',jl0 745 !WRITE(numout,*) ' vol = ',zvt_i_ini(ji,jj),SUM(zv_i_ini(:,ji,jj)) 746 ztest_1 = 0 747 ENDIF 748 749 ! Test 2: volume conservation 750 zV_cons = SUM(zv_i_ini(:,ji,jj)) 751 zconv = ABS(zvt_i_ini(ji,jj) - zV_cons) 752 753 IF ( zconv .LT. 1.0e-6 ) THEN 754 ztest_2 = 1 755 ELSE 756 ! this write is useful 757 !WRITE(numout,*) ' * TEST2 VOLUME NOT CONSERVED *** zV_cons = ', zV_cons, & 758 ! ' zvt_i_ini = ', zvt_i_ini(ji,jj) 759 !WRITE(numout,*) 'ji,jj,narea ',ji,jj,narea 760 !WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj) 761 !WRITE(numout,*) ' zhm_i_ini : ', zhm_i_ini(ji,jj) 762 !WRITE(numout,*) ' zht_i_ini(:,jij,jj) ',zht_i_ini(:,ji,jj) 763 !WRITE(numout,*) ' za_i_ini(:,jij,jj) ',za_i_ini(:,ji,jj) 764 !WRITE(numout,*) ' hi_max ',hi_max 765 !WRITE(numout,*) ' jl0 = ',jl0 766 ztest_2 = 0 767 ENDIF 768 769 ! Test 3: thickness of the last category is in-bounds ? 770 IF ( zht_i_ini(i_fill, ji,jj) .GT. hi_max(i_fill-1) ) THEN 771 ztest_3 = 1 772 ELSE 773 ! this write is useful 774 !WRITE(numout,*) ' * TEST 3 THICKNESS OF THE LAST CATEGORY OUT OF BOUNDS *** zht_i_ini(i_fill,ji,jj) = ', & 775 !zht_i_ini(i_fill,ji,jj), ' hi_max(jpl-1) = ', hi_max(i_fill-1) 776 !WRITE(numout,*) 'ji,jj,narea ',ji,jj,narea 777 !WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj) 778 !WRITE(numout,*) ' zhm_i_ini : ', zhm_i_ini(ji,jj) 779 !WRITE(numout,*) ' zht_i_ini(:,jij,jj) ',zht_i_ini(:,ji,jj) 780 !WRITE(numout,*) ' za_i_ini(:,jij,jj) ',za_i_ini(:,ji,jj) 781 !WRITE(numout,*) ' hi_max ',hi_max 782 !WRITE(numout,*) ' jl0 = ',jl0 783 ztest_3 = 0 784 ENDIF 785 786 ! Test 4: positivity of ice concentrations 787 ztest_4 = 1 788 DO jl = 1, jpl 789 IF ( za_i_ini(jl,ji,jj) .LT. 0._wp ) THEN 790 ! this write is useful 791 !WRITE(numout,*) ' * TEST 4 POSITIVITY NOT OK FOR CAT ', jl, 'WITH A = ', za_i_ini(jl,ji,jj) 792 !WRITE(numout,*) 'ji,jj,narea ',ji,jj,narea 793 !WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj) 794 !WRITE(numout,*) ' zhm_i_ini : ', zhm_i_ini(ji,jj) 795 !WRITE(numout,*) ' zht_i_ini(:,jij,jj) ',zht_i_ini(:,ji,jj) 796 !WRITE(numout,*) ' za_i_ini(:,jij,jj) ',za_i_ini(:,ji,jj) 797 !WRITE(numout,*) ' hi_max ',hi_max 798 !WRITE(numout,*) ' jl0 = ',jl0 799 !WRITE(numout,*) 800 ztest_4 = 0 801 ENDIF 802 END DO 803 804 ENDIF ! ztest_1 + ztest_2 + ztest_3 + ztest_4 805 806 ztests = ztest_1 + ztest_2 + ztest_3 + ztest_4 807 808 END DO ! i_fill 809 810 !WRITE(numout,*) ' ztests : ', ztests 811 !IF ( ztests .NE. 4 ) THEN 812 !WRITE(numout,*) 813 !WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ' 814 !WRITE(numout,*) ' !!!! RED ALERT !!! ' 815 !WRITE(numout,*) ' !!!! BIIIIP BIIIP BIIIIP BIIIIP !!!' 816 !WRITE(numout,*) ' !!!! Something is wrong in the LIM3 initialization procedure ' 817 !WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ' 818 !WRITE(numout,*) 'ji,jj,narea ',ji,jj,narea 819 !WRITE(numout,*) ' *** ztests is not equal to 4 ' 820 !WRITE(numout,*) ' *** ztest_i (i=1,4) = ', ztest_1, ztest_2,ztest_3,ztest_4 821 !WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj) 822 !WRITE(numout,*) ' zhm_i_ini : ', zhm_i_ini(ji,jj) 823 !WRITE(numout,*) ' zht_i_ini(:,jij,jj) ',zht_i_ini(:,ji,jj) 824 !WRITE(numout,*) ' za_i_ini(:,jij,jj) ',za_i_ini(:,ji,jj) 825 !WRITE(numout,*) ' hi_max ',hi_max 826 !ENDIF ! ztests .NE. 4 827 828 ENDIF ! jl0 ne 1 829 830 ENDIF ! zat_i_ini ne 0 831 END DO ! jj 832 END DO ! ji 833 834 835 !--------------------------------------------------------------------- 836 ! 3.3) Space-dependent arrays for ice state variables 837 !--------------------------------------------------------------------- 838 839 ! Ice concentration, thickness and volume, ice salinity, ice age, surface 840 ! temperature 841 DO jl = 1, jpl ! loop over categories 842 DO jj = 1, jpj 843 DO ji = 1, jpi 844 a_i(ji,jj,jl) = zswitch(ji,jj) * za_i_ini (jl,ji,jj) ! concentration 845 ht_i(ji,jj,jl) = zswitch(ji,jj) * zht_i_ini(jl,ji,jj) !ice thickness 846 847 IF( zhm_i_ini( ji,jj ) .GT. 0_wp )THEN ; ht_s(ji,jj,jl) = ht_i(ji,jj,jl) * ( zhm_s_ini( ji,jj ) / zhm_i_ini( ji,jj ) ) 848 ELSE ; ht_s(ji,jj,jl) = 0._wp 849 ENDIF 850 sm_i(ji,jj,jl) = zswitch(ji,jj) * zsm_i_ini(ji,jj) !+ (1._wp - zswitch(ji,jj) ) * rn_simin ! salinity 851 o_i(ji,jj,jl) = zswitch(ji,jj) * 1._wp + ( 1._wp -zswitch(ji,jj) ) ! age 852 t_su(ji,jj,jl) = sist_ini(ji,jj) 853 854 ! This case below should not be used if (ht_s/ht_i) is ok in 855 ! namelist 856 ! In case snow load is in excess that would lead to 857 ! transformation from snow to ice 858 ! Then, transfer the snow excess into the ice (different from 859 ! limthd_dh) 860 zdh = MAX( 0._wp, ( rhosn * ht_s(ji,jj,jl) + ( rhoic - rau0 ) *ht_i(ji,jj,jl) ) * r1_rau0 ) 861 ! recompute ht_i, ht_s avoiding out of bounds values 862 ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh ) 863 ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic *r1_rhosn ) 864 865 ! ice volume, salt content, age content 866 v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) !ice volume 867 v_s(ji,jj,jl) = ht_s(ji,jj,jl) * a_i(ji,jj,jl) !snow volume 868 smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) *v_i(ji,jj,jl) ! salt content 869 oa_i(ji,jj,jl) = o_i(ji,jj,jl) * a_i(ji,jj,jl) !age content 870 END DO ! ji 871 END DO ! jj 872 END DO ! jl 873 874 !cbr 875 DO jk = 1, nlay_s 876 DO jl = 1, jpl ! loop over categories 877 !rbb t_s(:,:,1,jl) = tbif_ini(:,:,1) 878 t_s(:,:,1,jl) = tbif_ini(:,:,1)*zswitch(:,:)+ ( 1._wp - zswitch(:,:) ) * rt0 879 END DO ! jl 880 END DO ! jk 881 882 ! Snow temperature and heat content 883 DO jk = 1, nlay_s 884 DO jl = 1, jpl ! loop over categories 885 DO jj = 1, jpj 886 DO ji = 1, jpi 887 !cbr??? t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 888 ! Snow energy of melting 889 e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus ) 890 891 ! Mutliply by volume, and divide by number of layers to get 892 ! heat content in J/m2 893 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) *r1_nlay_s 894 END DO ! ji 895 END DO ! jj 896 END DO ! jl 897 END DO ! jk 898 899 ! Ice salinity, temperature and heat content 900 DO jk = 1, nlay_i 901 DO jl = 1, jpl ! loop over categories 902 DO jj = 1, jpj 903 DO ji = 1, jpi 904 !cbr??? t_i(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 905 t_i(ji,jj,jk,jl) = tbif_ini(ji,jj,2)*zswitch(ji,jj)+ ( 1._wp - zswitch(ji,jj) ) * rt0 906 s_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(ji,jj) !+ ( 1._wp - zswitch(ji,jj) ) * rn_simin 907 ztmelts = - tmut * s_i(ji,jj,jk,jl) + rt0 !Melting temperature in K 908 909 ! heat content per unit volume 910 e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * ( cpic * ( ztmelts - t_i(ji,jj,jk,jl) ) & 911 + lfus * ( 1._wp - (ztmelts-rt0) /MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) & 912 - rcp * ( ztmelts - rt0 ) ) 913 914 ! Mutliply by ice volume, and divide by number of layers to 915 ! get heat content in J/m2 916 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i 917 END DO ! ji 918 END DO ! jj 919 END DO ! jl 920 END DO ! jk 921 922 !cbr tmp CALL wrk_dealloc(jpl,jpi,jpj, zht_i_ini, za_i_ini, zv_i_ini) 923 CALL wrk_dealloc(jpl,jpi,jpj, zv_i_ini) 924 CALL wrk_dealloc( jpl,jpi,jpj,zht_i_ini,za_i_ini) 925 CALL wrk_dealloc( jpi,jpj, zhm_i_ini, zat_i_ini, zvt_i_ini, zhm_s_ini,zsm_i_ini ) 926 CALL wrk_dealloc( jpi,jpj,zidto ) 927 928 END SUBROUTINE limini_file 929 930 931 INTEGER FUNCTION alloc_lim_istate_init() 932 !!----------------------------------------------------------------------------- 933 !! 934 !! 935 !! 936 !! 937 !!----------------------------------------------------------------------------- 938 INTEGER :: ierr(1) 939 !!----------------------------------------------------------------------------- 940 ALLOCATE( hicif_ini(jpi,jpj) , hsnif_ini(jpi,jpj) , frld_ini(jpi,jpj) , sist_ini(jpi,jpj) , zswitch(jpi,jpj) , tbif_ini(jpi,jpj,3) , Stat=ierr(1) ) 941 alloc_lim_istate_init = MAXVAL(ierr) 942 IF( alloc_lim_istate_init /= 0 ) CALL ctl_warn( 'lim_istate_init: failed to allocate arrays') 943 944 END FUNCTION alloc_lim_istate_init 513 945 #else 514 946 !!----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
