 Timestamp:
 20150325T10:01:32+01:00 (9 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
r3784 r5168 9 9 !! ! 200704 (A. Weaver) Merge with OPAVAR/NEMOVAR 10 10 !! NEMO 3.3 ! 201005 (D. Lea) Update to work with NEMO v3.2 11 !!  ! 201005 (D. Lea) add calc_month_len routine based on day_init 11 !!  ! 201005 (D. Lea) add calc_month_len routine based on day_init 12 12 !! 3.4 ! 201210 (A. Weaver and K. Mogensen) Fix for direct initialization 13 13 !! … … 43 43 IMPLICIT NONE 44 44 PRIVATE 45 45 46 PUBLIC asm_inc_rea_nam !: Read namelist for ASM 46 47 PUBLIC asm_inc_init !: Initialize the increment arrays and IAU weights 47 48 PUBLIC calc_date !: Compute the calendar date YYYYMMDD on a given step … … 70 71 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: u_bkg , v_bkg !: Background u & v velocity components 71 72 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: t_bkginc, s_bkginc !: Increment to the background T & S 72 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: u_bkginc, v_bkginc !: Increment to the u & vcomponents 73 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: u_bkginc, v_bkginc !: Increment to the u & vcomponents 73 74 REAL(wp), PUBLIC, DIMENSION(:) , ALLOCATABLE :: wgtiau !: IAU weights for each time step 74 75 #if defined key_asminc … … 78 79 INTEGER , PUBLIC :: nitbkg !: Time step of the background state used in the Jb term 79 80 INTEGER , PUBLIC :: nitdin !: Time step of the background state for direct initialization 80 INTEGER , PUBLIC :: nitiaustr !: Time step of the start of the IAU interval 81 INTEGER , PUBLIC :: nitiaustr !: Time step of the start of the IAU interval 81 82 INTEGER , PUBLIC :: nitiaufin !: Time step of the end of the IAU interval 82 ! 83 ! 83 84 INTEGER , PUBLIC :: niaufn !: Type of IAU weighing function: = 0 Constant weighting 84 ! !: = 1 Linear hatlike, centred in middle of IAU interval 85 ! !: = 1 Linear hatlike, centred in middle of IAU interval 85 86 REAL(wp), PUBLIC :: salfixmin !: Ensure that the salinity is larger than this value if (ln_salfix) 86 87 … … 100 101 CONTAINS 101 102 102 SUBROUTINE asm_inc_init 103 !! 104 !! *** ROUTINE asm_inc_init *** 105 !! 106 !! ** Purpose : Initialize the assimilation increment and IAU weights. 107 !! 108 !! ** Method : Initialize the assimilation increment and IAU weights. 109 !! 110 !! ** Action : 111 !! 112 !! 113 !! 114 INTEGER :: ji,jj,jk 115 INTEGER :: jt 116 INTEGER :: imid 117 INTEGER :: inum 118 INTEGER :: iiauper ! Number of time steps in the IAU period 119 INTEGER :: icycper ! Number of time steps in the cycle 120 INTEGER :: iitend_date ! Date YYYYMMDD of final time step 121 INTEGER :: iitbkg_date ! Date YYYYMMDD of background time step for Jb term 122 INTEGER :: iitdin_date ! Date YYYYMMDD of background time step for DI 123 INTEGER :: iitiaustr_date ! Date YYYYMMDD of IAU interval start time step 124 INTEGER :: iitiaufin_date ! Date YYYYMMDD of IAU interval final time step 125 126 REAL(wp) :: znorm ! Normalization factor for IAU weights 127 REAL(wp) :: ztotwgt ! Value of timeintegrated IAU weights 128 ! (should be equal to one) 129 REAL(wp) :: z_inc_dateb ! Start date of interval on which increment is valid 130 REAL(wp) :: z_inc_datef ! End date of interval on which increment is valid 131 REAL(wp) :: zdate_bkg ! Date in background state file for DI 132 REAL(wp) :: zdate_inc ! Time axis in increments file 133 134 REAL(wp), POINTER, DIMENSION(:,:) :: hdiv 135 !! 103 SUBROUTINE asm_inc_rea_nam 136 104 NAMELIST/nam_asminc/ ln_bkgwri, & 137 105 & ln_trainc, ln_dyninc, ln_sshinc, & … … 158 126 salfixmin = 9999 159 127 nitbkg = 0 160 nitdin = 0 128 nitdin = 0 161 129 nitiaustr = 1 162 130 nitiaufin = 150 ! = 10 days with ORCA2 … … 187 155 WRITE(numout,*) ' Minimum salinity after applying the increments salfixmin = ', salfixmin 188 156 ENDIF 157 END SUBROUTINE 158 159 SUBROUTINE asm_inc_init 160 !! 161 !! *** ROUTINE asm_inc_init *** 162 !! 163 !! ** Purpose : Initialize the assimilation increment and IAU weights. 164 !! 165 !! ** Method : Initialize the assimilation increment and IAU weights. 166 !! 167 !! ** Action : 168 !! 169 !! 170 !! 171 INTEGER :: ji,jj,jk 172 INTEGER :: jt 173 INTEGER :: imid 174 INTEGER :: inum 175 INTEGER :: iiauper ! Number of time steps in the IAU period 176 INTEGER :: icycper ! Number of time steps in the cycle 177 INTEGER :: iitend_date ! Date YYYYMMDD of final time step 178 INTEGER :: iitbkg_date ! Date YYYYMMDD of background time step for Jb term 179 INTEGER :: iitdin_date ! Date YYYYMMDD of background time step for DI 180 INTEGER :: iitiaustr_date ! Date YYYYMMDD of IAU interval start time step 181 INTEGER :: iitiaufin_date ! Date YYYYMMDD of IAU interval final time step 182 183 REAL(wp) :: znorm ! Normalization factor for IAU weights 184 REAL(wp) :: ztotwgt ! Value of timeintegrated IAU weights 185 ! (should be equal to one) 186 REAL(wp) :: z_inc_dateb ! Start date of interval on which increment is valid 187 REAL(wp) :: z_inc_datef ! End date of interval on which increment is valid 188 REAL(wp) :: zdate_bkg ! Date in background state file for DI 189 REAL(wp) :: zdate_inc ! Time axis in increments file 190 191 REAL(wp), POINTER, DIMENSION(:,:) :: hdiv 192 !! 189 193 190 194 nitbkg_r = nitbkg + nit000  1 ! Background time referenced to nit000 … … 293 297 294 298 ! 295 ! Constant IAU forcing 299 ! Constant IAU forcing 296 300 ! 297 301 … … 303 307 304 308 ! 305 ! Linear hatlike, centred in middle of IAU interval 309 ! Linear hatlike, centred in middle of IAU interval 306 310 ! 307 311 … … 309 313 znorm = 0.0 310 314 IF ( MOD( iiauper, 2 ) == 0 ) THEN ! Even number of time steps in IAU interval 311 imid = iiauper / 2 315 imid = iiauper / 2 312 316 DO jt = 1, imid 313 317 znorm = znorm + REAL( jt ) … … 315 319 znorm = 2.0 * znorm 316 320 ELSE ! Odd number of time steps in IAU interval 317 imid = ( iiauper + 1 ) / 2 321 imid = ( iiauper + 1 ) / 2 318 322 DO jt = 1, imid  1 319 323 znorm = znorm + REAL( jt ) … … 342 346 DO jt = 1, icycper 343 347 ztotwgt = ztotwgt + wgtiau(jt) 344 WRITE(numout,*) ' ', jt, ' ', wgtiau(jt) 345 END DO 348 WRITE(numout,*) ' ', jt, ' ', wgtiau(jt) 349 END DO 346 350 WRITE(numout,*) ' ===================================' 347 351 WRITE(numout,*) ' Timeintegrated weight = ', ztotwgt 348 352 WRITE(numout,*) ' ===================================' 349 353 ENDIF 350 354 351 355 ENDIF 352 356 … … 381 385 CALL iom_open( c_asminc, inum ) 382 386 383 CALL iom_get( inum, 'time', zdate_inc ) 387 CALL iom_get( inum, 'time', zdate_inc ) 384 388 385 389 CALL iom_get( inum, 'z_inc_dateb', z_inc_dateb ) … … 389 393 390 394 IF(lwp) THEN 391 WRITE(numout,*) 395 WRITE(numout,*) 392 396 WRITE(numout,*) 'asm_inc_init : Assimilation increments valid ', & 393 397 & ' between dates ', NINT( z_inc_dateb ),' and ', & … … 405 409 & ' not agree with Direct Initialization time' ) 406 410 407 IF ( ln_trainc ) THEN 411 IF ( ln_trainc ) THEN 408 412 CALL iom_get( inum, jpdom_autoglo, 'bckint', t_bkginc, 1 ) 409 413 CALL iom_get( inum, jpdom_autoglo, 'bckins', s_bkginc, 1 ) … … 417 421 ENDIF 418 422 419 IF ( ln_dyninc ) THEN 420 CALL iom_get( inum, jpdom_autoglo, 'bckinu', u_bkginc, 1 ) 421 CALL iom_get( inum, jpdom_autoglo, 'bckinv', v_bkginc, 1 ) 423 IF ( ln_dyninc ) THEN 424 CALL iom_get( inum, jpdom_autoglo, 'bckinu', u_bkginc, 1 ) 425 CALL iom_get( inum, jpdom_autoglo, 'bckinv', v_bkginc, 1 ) 422 426 ! Apply the masks 423 427 u_bkginc(:,:,:) = u_bkginc(:,:,:) * umask(:,:,:) … … 428 432 WHERE( ABS( v_bkginc(:,:,:) ) > 1.0e+10 ) v_bkginc(:,:,:) = 0.0 429 433 ENDIF 430 434 431 435 IF ( ln_sshinc ) THEN 432 436 CALL iom_get( inum, jpdom_autoglo, 'bckineta', ssh_bkginc, 1 ) … … 448 452 449 453 CALL iom_close( inum ) 450 454 451 455 ENDIF 452 456 … … 457 461 IF ( ln_dyninc .AND. nn_divdmp > 0 ) THEN 458 462 459 CALL wrk_alloc(jpi,jpj,hdiv) 463 CALL wrk_alloc(jpi,jpj,hdiv) 460 464 461 465 DO jt = 1, nn_divdmp … … 482 486 u_bkginc(ji,jj,jk) = u_bkginc(ji,jj,jk) + 0.2_wp * ( e1t(ji+1,jj)*e2t(ji+1,jj) * hdiv(ji+1,jj) & 483 487  e1t(ji ,jj)*e2t(ji ,jj) * hdiv(ji ,jj) ) & 484 / e1u(ji,jj) * umask(ji,jj,jk) 488 / e1u(ji,jj) * umask(ji,jj,jk) 485 489 v_bkginc(ji,jj,jk) = v_bkginc(ji,jj,jk) + 0.2_wp * ( e1t(ji,jj+1)*e2t(ji,jj+1) * hdiv(ji,jj+1) & 486 490  e1t(ji,jj )*e2t(ji,jj ) * hdiv(ji,jj ) ) & 487 / e2v(ji,jj) * vmask(ji,jj,jk) 491 / e2v(ji,jj) * vmask(ji,jj,jk) 488 492 END DO 489 493 END DO … … 493 497 END DO 494 498 495 CALL wrk_dealloc(jpi,jpj,hdiv) 499 CALL wrk_dealloc(jpi,jpj,hdiv) 496 500 497 501 ENDIF … … 523 527 CALL iom_open( c_asmdin, inum ) 524 528 525 CALL iom_get( inum, 'rdastp', zdate_bkg ) 526 529 CALL iom_get( inum, 'rdastp', zdate_bkg ) 530 527 531 IF(lwp) THEN 528 WRITE(numout,*) 532 WRITE(numout,*) 529 533 WRITE(numout,*) 'asm_inc_init : Assimilation background state valid at : ', & 530 534 & NINT( zdate_bkg ) … … 536 540 & ' not agree with Direct Initialization time' ) 537 541 538 IF ( ln_trainc ) THEN 542 IF ( ln_trainc ) THEN 539 543 CALL iom_get( inum, jpdom_autoglo, 'tn', t_bkg ) 540 544 CALL iom_get( inum, jpdom_autoglo, 'sn', s_bkg ) … … 543 547 ENDIF 544 548 545 IF ( ln_dyninc ) THEN 549 IF ( ln_dyninc ) THEN 546 550 CALL iom_get( inum, jpdom_autoglo, 'un', u_bkg ) 547 551 CALL iom_get( inum, jpdom_autoglo, 'vn', v_bkg ) … … 549 553 v_bkg(:,:,:) = v_bkg(:,:,:) * vmask(:,:,:) 550 554 ENDIF 551 555 552 556 IF ( ln_sshinc ) THEN 553 557 CALL iom_get( inum, jpdom_autoglo, 'sshn', ssh_bkg ) … … 565 569 !! 566 570 !! *** ROUTINE calc_date *** 567 !! 571 !! 568 572 !! ** Purpose : Compute the calendar date YYYYMMDD at a given time step. 569 573 !! 570 574 !! ** Method : Compute the calendar date YYYYMMDD at a given time step. 571 575 !! 572 !! ** Action : 576 !! ** Action : 573 577 !! 574 578 INTEGER, INTENT(IN) :: kit000 ! Initial time step … … 595 599 iyea0 = kdate0 / 10000 596 600 imon0 = ( kdate0  ( iyea0 * 10000 ) ) / 100 597 iday0 = kdate0  ( iyea0 * 10000 )  ( imon0 * 100 ) 601 iday0 = kdate0  ( iyea0 * 10000 )  ( imon0 * 100 ) 598 602 599 603 ! Check that kt >= kit000  1 … … 610 614 ! Compute the number of days from the initial date 611 615 idaystp = INT( REAL( kt  kit000 ) * rdt / 86400. ) 612 616 613 617 iday = iday0 614 618 imon = imon0 … … 628 632 iyea = iyea + 1 629 633 CALL calc_month_len( iyea, imonth_len ) ! update month lengths 630 ENDIF 634 ENDIF 631 635 idaycnt = idaycnt + 1 632 636 END DO … … 640 644 !! 641 645 !! *** ROUTINE calc_month_len *** 642 !! 646 !! 643 647 !! ** Purpose : Compute the number of days in a months given a year. 644 648 !! 645 !! ** Method : 649 !! ** Method : 646 650 !! 647 651 INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year … … 650 654 ! 651 655 ! length of the month of the current year (from nleapy, read in namelist) 652 IF ( nleapy < 2 ) THEN 656 IF ( nleapy < 2 ) THEN 653 657 imonth_len(:) = (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /) 654 658 IF ( nleapy == 1 ) THEN ! we are using calendar with leap years … … 667 671 !! 668 672 !! *** ROUTINE tra_asm_inc *** 669 !! 673 !! 670 674 !! ** Purpose : Apply the tracer (T and S) assimilation increments 671 675 !! 672 676 !! ** Method : Direct initialization or Incremental Analysis Updating 673 677 !! 674 !! ** Action : 678 !! ** Action : 675 679 !! 676 680 INTEGER, INTENT(IN) :: kt ! Current time step … … 682 686 !! 683 687 684 ! freezing point calculation taken from oc_fz_pt (but calculated for all depths) 685 ! used to prevent the applied increments taking the temperature below the local freezing point 686 687 #if defined key_cice 688 ! freezing point calculation taken from oc_fz_pt (but calculated for all depths) 689 ! used to prevent the applied increments taking the temperature below the local freezing point 690 691 #if defined key_cice 688 692 fzptnz(:,:,:) = 1.8_wp 689 #else 693 #else 690 694 DO jk = 1, jpk 691 695 DO jj = 1, jpj 692 696 DO ji = 1, jpk 693 fzptnz (ji,jj,jk) = ( 0.0575_wp + 1.710523e3_wp * SQRT( tsn(ji,jj,jk,jp_sal) ) & 694  2.154996e4_wp * tsn(ji,jj,jk,jp_sal) ) * tsn(ji,jj,jk,jp_sal) & 695  7.53e4_wp * fsdepw(ji,jj,jk) ! (pressure in dbar) 697 fzptnz (ji,jj,jk) = ( 0.0575_wp + 1.710523e3_wp * SQRT( tsn(ji,jj,jk,jp_sal) ) & 698  2.154996e4_wp * tsn(ji,jj,jk,jp_sal) ) * tsn(ji,jj,jk,jp_sal) & 699  7.53e4_wp * fsdepw(ji,jj,jk) ! (pressure in dbar) 696 700 END DO 697 701 END DO 698 702 END DO 699 #endif 703 #endif 700 704 701 705 IF ( ln_asmiau ) THEN … … 711 715 712 716 IF(lwp) THEN 713 WRITE(numout,*) 717 WRITE(numout,*) 714 718 WRITE(numout,*) 'tra_asm_inc : Tracer IAU at time step = ', & 715 719 & kt,' with IAU weight = ', wgtiau(it) … … 722 726 ! Do not apply negative increments if the temperature will fall below freezing 723 727 WHERE(t_bkginc(:,:,jk) > 0.0_wp .OR. & 724 & t sn(:,:,jk,jp_tem) + tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * wgtiau(it) > fzptnz(:,:,jk) )725 tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt 728 & t_bkg(:,:,jk) + tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * wgtiau(it) > fzptnz(:,:,jk) ) 729 tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt 726 730 END WHERE 727 731 ELSE 728 tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt 732 tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt 729 733 ENDIF 730 734 IF (ln_salfix) THEN … … 732 736 ! minimum value salfixmin 733 737 WHERE(s_bkginc(:,:,jk) > 0.0_wp .OR. & 734 & tsn(:,:,jk,jp_sal) + tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * wgtiau(it) > salfixmin )738 & s_bkg(:,:,jk) + tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * wgtiau(it) > salfixmin ) 735 739 tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt 736 740 END WHERE … … 753 757 ! Direct Initialization 754 758 ! 755 759 756 760 IF ( kt == nitdin_r ) THEN 757 761 … … 762 766 ! Do not apply negative increments if the temperature will fall below freezing 763 767 WHERE(t_bkginc(:,:,:) > 0.0_wp .OR. & 764 & tsn(:,:,:,jp_tem) + t_bkginc(:,:,:) > fzptnz(:,:,:) ) 765 tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:) 768 & tsn(:,:,:,jp_tem) + t_bkginc(:,:,:) > fzptnz(:,:,:) ) 769 tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:) 766 770 END WHERE 767 771 ELSE 768 tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:) 772 tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:) 769 773 ENDIF 770 774 IF (ln_salfix) THEN … … 772 776 ! minimum value salfixmin 773 777 WHERE(s_bkginc(:,:,:) > 0.0_wp .OR. & 774 & tsn(:,:,:,jp_sal) + s_bkginc(:,:,:) > salfixmin ) 775 tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:) 778 & tsn(:,:,:,jp_sal) + s_bkginc(:,:,:) > salfixmin ) 779 tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:) 776 780 END WHERE 777 781 ELSE 778 tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:) 782 tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:) 779 783 ENDIF 780 784 … … 782 786 783 787 CALL eos( tsb, rhd, rhop ) ! Before potential and in situ densities 784 788 785 789 IF( ln_zps .AND. .NOT. lk_c1d ) & 786 790 & CALL zps_hde( nit000, jpts, tsb, & ! Partial steps: before horizontal derivative … … 797 801 DEALLOCATE( s_bkg ) 798 802 ENDIF 799 ! 803 ! 800 804 ENDIF 801 805 ! Perhaps the following call should be in step … … 808 812 !! 809 813 !! *** ROUTINE dyn_asm_inc *** 810 !! 814 !! 811 815 !! ** Purpose : Apply the dynamics (u and v) assimilation increments. 812 816 !! 813 817 !! ** Method : Direct initialization or Incremental Analysis Updating. 814 818 !! 815 !! ** Action : 819 !! ** Action : 816 820 !! 817 821 INTEGER, INTENT(IN) :: kt ! Current time step … … 834 838 835 839 IF(lwp) THEN 836 WRITE(numout,*) 840 WRITE(numout,*) 837 841 WRITE(numout,*) 'dyn_asm_inc : Dynamics IAU at time step = ', & 838 842 & kt,' with IAU weight = ', wgtiau(it) … … 845 849 va(:,:,jk) = va(:,:,jk) + v_bkginc(:,:,jk) * zincwgt 846 850 END DO 847 851 848 852 IF ( kt == nitiaufin_r ) THEN 849 853 DEALLOCATE( u_bkginc ) … … 853 857 ENDIF 854 858 855 ELSEIF ( ln_asmdin ) THEN 859 ELSEIF ( ln_asmdin ) THEN 856 860 857 861 ! 858 862 ! Direct Initialization 859 863 ! 860 864 861 865 IF ( kt == nitdin_r ) THEN 862 866 … … 865 869 ! Initialize the now fields with the background + increment 866 870 un(:,:,:) = u_bkg(:,:,:) + u_bkginc(:,:,:) 867 vn(:,:,:) = v_bkg(:,:,:) + v_bkginc(:,:,:) 871 vn(:,:,:) = v_bkg(:,:,:) + v_bkginc(:,:,:) 868 872 869 873 ub(:,:,:) = un(:,:,:) ! Update before fields 870 874 vb(:,:,:) = vn(:,:,:) 871 875 872 876 DEALLOCATE( u_bkg ) 873 877 DEALLOCATE( v_bkg ) … … 884 888 !! 885 889 !! *** ROUTINE ssh_asm_inc *** 886 !! 890 !! 887 891 !! ** Purpose : Apply the sea surface height assimilation increment. 888 892 !! 889 893 !! ** Method : Direct initialization or Incremental Analysis Updating. 890 894 !! 891 !! ** Action : 895 !! ** Action : 892 896 !! 893 897 INTEGER, INTENT(IN) :: kt ! Current time step … … 910 914 911 915 IF(lwp) THEN 912 WRITE(numout,*) 916 WRITE(numout,*) 913 917 WRITE(numout,*) 'ssh_asm_inc : SSH IAU at time step = ', & 914 918 & kt,' with IAU weight = ', wgtiau(it) … … 938 942 939 943 ! Initialize the now fields the background + increment 940 sshn(:,:) = ssh_bkg(:,:) + ssh_bkginc(:,:) 944 sshn(:,:) = ssh_bkg(:,:) + ssh_bkginc(:,:) 941 945 942 946 ! Update before fields 943 sshb(:,:) = sshn(:,:) 947 sshb(:,:) = sshn(:,:) 944 948 945 949 IF( lk_vvl ) THEN … … 961 965 !! 962 966 !! *** ROUTINE seaice_asm_inc *** 963 !! 967 !! 964 968 !! ** Purpose : Apply the sea ice assimilation increment. 965 969 !! 966 970 !! ** Method : Direct initialization or Incremental Analysis Updating. 967 971 !! 968 !! ** Action : 972 !! ** Action : 969 973 !! 970 974 !! History : … … 998 1002 999 1003 it = kt  nit000 + 1 1000 zincwgt = wgtiau(it) ! IAU weight for the current time step 1004 zincwgt = wgtiau(it) ! IAU weight for the current time step 1001 1005 ! note this is not a tendency so should not be divided by rdt (as with the tracer and other increments) 1002 1006 1003 1007 IF(lwp) THEN 1004 WRITE(numout,*) 1008 WRITE(numout,*) 1005 1009 WRITE(numout,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', & 1006 1010 & kt,' with IAU weight = ', wgtiau(it) … … 1021 1025 ! Nudge sea ice depth to bring it up to a required minimum depth 1022 1026 1023 WHERE( zseaicendg(:,:) > 0.0_wp .AND. hicif(:,:) < zhicifmin ) 1024 zhicifinc(:,:) = (zhicifmin  hicif(:,:)) * zincwgt 1027 WHERE( zseaicendg(:,:) > 0.0_wp .AND. hicif(:,:) < zhicifmin ) 1028 zhicifinc(:,:) = (zhicifmin  hicif(:,:)) * zincwgt 1025 1029 ELSEWHERE 1026 1030 zhicifinc(:,:) = 0.0_wp … … 1029 1033 ! nudge ice depth 1030 1034 hicif(:,:)=hicif(:,:) + zhicifinc(:,:) 1031 phicif(:,:)=phicif(:,:) + zhicifinc(:,:) 1035 phicif(:,:)=phicif(:,:) + zhicifinc(:,:) 1032 1036 1033 1037 ! seaice salinity balancing (to add) … … 1071 1075 zofrld(:,:)=frld(:,:) 1072 1076 zohicif(:,:)=hicif(:,:) 1073 1077 1074 1078 ! Initialize the now fields the background + increment 1075 1079 1076 1080 frld(:,:) = MIN( MAX( frld(:,:)  seaice_bkginc(:,:), 0.0_wp), 1.0_wp) 1077 pfrld(:,:) = frld(:,:) 1081 pfrld(:,:) = frld(:,:) 1078 1082 fr_i(:,:) = 1.0_wp  frld(:,:) ! adjust ice fraction 1079 1083 … … 1082 1086 ! Nudge sea ice depth to bring it up to a required minimum depth 1083 1087 1084 WHERE( zseaicendg(:,:) > 0.0_wp .AND. hicif(:,:) < zhicifmin ) 1085 zhicifinc(:,:) = (zhicifmin  hicif(:,:)) * zincwgt 1088 WHERE( zseaicendg(:,:) > 0.0_wp .AND. hicif(:,:) < zhicifmin ) 1089 zhicifinc(:,:) = (zhicifmin  hicif(:,:)) * zincwgt 1086 1090 ELSEWHERE 1087 1091 zhicifinc(:,:) = 0.0_wp … … 1090 1094 ! nudge ice depth 1091 1095 hicif(:,:)=hicif(:,:) + zhicifinc(:,:) 1092 phicif(:,:)=phicif(:,:) 1096 phicif(:,:)=phicif(:,:) 1093 1097 1094 1098 ! seaice salinity balancing (to add) 1095 1096 #endif 1097 1099 1100 #endif 1101 1098 1102 #if defined key_cice 1099 1103 … … 1110 1114 #if defined key_cice 1111 1115 1112 ! Zero ice increment tendency into CICE 1116 ! Zero ice increment tendency into CICE 1113 1117 ndaice_da(:,:) = 0.0_wp 1114 1118 1115 1119 #endif 1116 1120 1117 1121 ENDIF 1118 1122 1119 1123 !#if defined key_lim2  defined key_cice 1120 1124 ! 1121 ! IF (ln_seaicebal ) THEN 1125 ! IF (ln_seaicebal ) THEN 1122 1126 ! !! balancing salinity increments 1123 1127 ! !! simple case from limflx.F90 (doesn't include a mass flux) … … 1131 1135 ! 1132 1136 ! DO jj = 1, jpj 1133 ! DO ji = 1, jpi 1137 ! DO ji = 1, jpi 1134 1138 ! ! calculate change in ice and snow mass per unit area 1135 1139 ! ! positive values imply adding salt to the ocean (results from ice formation) … … 1142 1146 ! 1143 1147 ! ! prevent small mld 1144 ! ! less than 10m can cause salinity instability 1148 ! ! less than 10m can cause salinity instability 1145 1149 ! IF (mld < 10) mld=10 1146 1150 ! 1147 ! ! set to bottom of a level 1151 ! ! set to bottom of a level 1148 1152 ! DO jk = jpk1, 2, 1 1149 ! IF ((mld > gdepw(ji,jj,jk)) .and. (mld < gdepw(ji,jj,jk+1))) THEN 1153 ! IF ((mld > gdepw(ji,jj,jk)) .and. (mld < gdepw(ji,jj,jk+1))) THEN 1150 1154 ! mld=gdepw(ji,jj,jk+1) 1151 1155 ! jkmax=jk … … 1154 1158 ! 1155 1159 ! ! avoid applying salinity balancing in shallow water or on land 1156 ! ! 1160 ! ! 1157 1161 ! 1158 1162 ! ! dsal_ocn (psu kg m^2) / (kg m^3 * m) … … 1165 1169 ! 1166 1170 ! ! put increments in for levels in the mixed layer 1167 ! ! but prevent salinity below a threshold value 1168 ! 1169 ! DO jk = 1, jkmax 1170 ! 1171 ! IF (dsal_ocn > 0.0_wp .or. sb(ji,jj,jk)+dsal_ocn > sal_thresh) THEN 1171 ! ! but prevent salinity below a threshold value 1172 ! 1173 ! DO jk = 1, jkmax 1174 ! 1175 ! IF (dsal_ocn > 0.0_wp .or. sb(ji,jj,jk)+dsal_ocn > sal_thresh) THEN 1172 1176 ! sb(ji,jj,jk) = sb(ji,jj,jk) + dsal_ocn 1173 1177 ! sn(ji,jj,jk) = sn(ji,jj,jk) + dsal_ocn … … 1180 1184 ! ! 1181 1185 ! !! Adjust fsalt. A +ve fsalt means adding salt to ocean 1182 ! !! fsalt(ji,jj) = fsalt(ji,jj) + zpmess ! adjust fsalt 1183 ! !! 1184 ! !! emps(ji,jj) = emps(ji,jj) + zpmess ! or adjust emps (see icestp1d) 1186 ! !! fsalt(ji,jj) = fsalt(ji,jj) + zpmess ! adjust fsalt 1187 ! !! 1188 ! !! emps(ji,jj) = emps(ji,jj) + zpmess ! or adjust emps (see icestp1d) 1185 1189 ! !! ! EP (kg m2 s2) 1186 1190 ! ! emp(ji,jj) = emp(ji,jj) + zpmess ! EP (kg m2 s2)
Note: See TracChangeset
for help on using the changeset viewer.