Changeset 5168
- Timestamp:
- 2015-03-25T10:01:32+01:00 (8 years ago)
- Location:
- branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO
- Files:
-
- 4 edited
-
OPATAM_SRC/tstool_tam.F90 (modified) (1 diff)
-
OPA_SRC/ASM/asminc.F90 (modified) (62 diffs)
-
OPA_SRC/DOM/istate.F90 (modified) (28 diffs)
-
OPA_SRC/nemogcm.F90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPATAM_SRC/tstool_tam.F90
r3611 r5168 14 14 & e1u, & 15 15 & e2u, & 16 #if defined key_zco17 16 & e3t_0, & 18 #else19 17 & e3u, & 20 #endif21 18 & umask, & 22 19 & mig, & -
branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
r3784 r5168 9 9 !! ! 2007-04 (A. Weaver) Merge with OPAVAR/NEMOVAR 10 10 !! NEMO 3.3 ! 2010-05 (D. Lea) Update to work with NEMO v3.2 11 !! - ! 2010-05 (D. Lea) add calc_month_len routine based on day_init 11 !! - ! 2010-05 (D. Lea) add calc_month_len routine based on day_init 12 12 !! 3.4 ! 2012-10 (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- & v-components 73 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: u_bkginc, v_bkginc !: Increment to the u- & v-components 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 hat-like, centred in middle of IAU interval 85 ! !: = 1 Linear hat-like, 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 time-integrated 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 time-integrated 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 hat-like, centred in middle of IAU interval 309 ! Linear hat-like, 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,*) ' Time-integrated 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.710523e-3_wp * SQRT( tsn(ji,jj,jk,jp_sal) ) & 694 - 2.154996e-4_wp * tsn(ji,jj,jk,jp_sal) ) * tsn(ji,jj,jk,jp_sal) & 695 - 7.53e-4_wp * fsdepw(ji,jj,jk) ! (pressure in dbar) 697 fzptnz (ji,jj,jk) = ( -0.0575_wp + 1.710523e-3_wp * SQRT( tsn(ji,jj,jk,jp_sal) ) & 698 - 2.154996e-4_wp * tsn(ji,jj,jk,jp_sal) ) * tsn(ji,jj,jk,jp_sal) & 699 - 7.53e-4_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 = jpk-1, 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 ! !! ! E-P (kg m-2 s-2) 1186 1190 ! ! emp(ji,jj) = emp(ji,jj) + zpmess ! E-P (kg m-2 s-2) -
branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
r3711 r5168 10 10 !! 8.0 ! 2001-09 (M. Levy, M. Ben Jelloul) istate_uvg 11 11 !! NEMO 1.0 ! 2003-08 (G. Madec, C. Talandier) F90: Free form, modules + EEL R5 12 !! - ! 2004-05 (A. Koch-Larrouy) istate_gyre 12 !! - ! 2004-05 (A. Koch-Larrouy) istate_gyre 13 13 !! 2.0 ! 2006-07 (S. Masson) distributed restart using iom 14 14 !! 3.3 ! 2010-10 (C. Ethe) merge TRC-TRA 15 !! 3.4 ! 2011-04 (G. Madec) Merge of dtatem and dtasal & suppression of tb,tn/sb,sn 15 !! 3.4 ! 2011-04 (G. Madec) Merge of dtatem and dtasal & suppression of tb,tn/sb,sn 16 16 !!---------------------------------------------------------------------- 17 17 … … 24 24 !! istate_uvg : initial velocity in geostropic balance 25 25 !!---------------------------------------------------------------------- 26 USE oce ! ocean dynamics and active tracers 27 USE dom_oce ! ocean space and time domain 26 USE oce ! ocean dynamics and active tracers 27 USE dom_oce ! ocean space and time domain 28 28 USE daymod ! calendar 29 29 USE eosbn2 ! eq. of state, Brunt Vaisala frequency (eos routine) … … 46 46 USE wrk_nemo ! Memory allocation 47 47 USE timing ! Timing 48 USE asminc 48 49 49 50 IMPLICIT NONE … … 65 66 !!---------------------------------------------------------------------- 66 67 !! *** ROUTINE istate_init *** 67 !! 68 !! 68 69 !! ** Purpose : Initialization of the dynamics and tracer fields. 69 70 !!---------------------------------------------------------------------- … … 83 84 rhd (:,:,: ) = 0.e0 84 85 rhop (:,:,: ) = 0.e0 85 rn2 (:,:,: ) = 0.e0 86 tsa (:,:,:,:) = 0.e0 86 rn2 (:,:,: ) = 0.e0 87 tsa (:,:,:,:) = 0.e0 87 88 88 89 IF( ln_rstart ) THEN ! Restart from a file … … 100 101 CALL day_init ! model calendar (using both namelist and restart infos) 101 102 ! ! Initialization of ocean to zero 102 ! before fields ! now fields 103 ! before fields ! now fields 103 104 sshb (:,:) = 0._wp ; sshn (:,:) = 0._wp 104 105 ub (:,:,:) = 0._wp ; un (:,:,:) = 0._wp 105 vb (:,:,:) = 0._wp ; vn (:,:,:) = 0._wp 106 vb (:,:,:) = 0._wp ; vn (:,:,:) = 0._wp 106 107 rotb (:,:,:) = 0._wp ; rotn (:,:,:) = 0._wp 107 108 hdivb(:,:,:) = 0._wp ; hdivn(:,:,:) = 0._wp … … 109 110 IF( cp_cfg == 'eel' ) THEN 110 111 CALL istate_eel ! EEL configuration : start from pre-defined U,V T-S fields 111 ELSEIF( cp_cfg == 'gyre' ) THEN 112 ELSEIF( cp_cfg == 'gyre' ) THEN 112 113 CALL istate_gyre ! GYRE configuration : start from pre-defined T-S fields 113 114 ELSEIF( ln_tsd_init ) THEN ! Initial T-S fields read in files … … 124 125 & rhd, gru , grv ) ! of t,s,rd at ocean bottom 125 126 #endif 126 ! 127 ! 127 128 ! - ML - sshn could be modified by istate_eel, so that initialization of fse3t_b is done here 128 129 IF( lk_vvl ) THEN … … 133 134 ! ! define e3u_b, e3v_b from e3t_b initialized in domzgr 134 135 CALL dom_vvl_2( nit000, fse3u_b(:,:,:), fse3v_b(:,:,:) ) 135 ! 136 ! 136 137 ENDIF 137 138 ! … … 146 147 ENDIF 147 148 ! 149 IF ( lk_asminc .AND. ln_asmdin ) THEN 150 neuler = 0 ! restart with an euler from the corrected background 151 ENDIF 152 148 153 IF( nn_timing == 1 ) CALL timing_stop('istate_init') 149 154 ! … … 153 158 !!--------------------------------------------------------------------- 154 159 !! *** ROUTINE istate_t_s *** 155 !! 156 !! ** Purpose : Intialization of the temperature field with an 160 !! 161 !! ** Purpose : Intialization of the temperature field with an 157 162 !! analytical profile or a file (i.e. in EEL configuration) 158 163 !! … … 184 189 !!---------------------------------------------------------------------- 185 190 !! *** ROUTINE istate_eel *** 186 !! 191 !! 187 192 !! ** Purpose : Initialization of the dynamics and tracers for EEL R5 188 193 !! configuration (channel with or without a topographic bump) … … 195 200 USE divcur ! hor. divergence & rel. vorticity (div_cur routine) 196 201 USE iom 197 202 198 203 INTEGER :: inum ! temporary logical unit 199 204 INTEGER :: ji, jj, jk ! dummy loop indices … … 207 212 !!---------------------------------------------------------------------- 208 213 209 SELECT CASE ( jp_cfg ) 214 SELECT CASE ( jp_cfg ) 210 215 ! ! ==================== 211 216 CASE ( 5 ) ! EEL R5 configuration … … 244 249 ! set the dynamics: U,V, hdiv, rot (and ssh if necessary) 245 250 ! ---------------- 246 ! Start EEL5 configuration with barotropic geostrophic velocities 251 ! Start EEL5 configuration with barotropic geostrophic velocities 247 252 ! according the sshb and sshn SSH imposed. 248 253 ! we assume a uniform grid (hence the use of e1t(1,1) for delta_y) 249 ! we use the Coriolis frequency at mid-channel. 254 ! we use the Coriolis frequency at mid-channel. 250 255 ub(:,:,:) = zueel * umask(:,:,:) 251 256 un(:,:,:) = ub(:,:,:) … … 254 259 ! 255 260 DO jj = 1, jpjglo 256 zssh(:,jj) = - (FLOAT(jj)- FLOAT(jpjglo-1)/2.)*zueel*e1t(1,1)*zfcor/grav 261 zssh(:,jj) = - (FLOAT(jj)- FLOAT(jpjglo-1)/2.)*zueel*e1t(1,1)*zfcor/grav 257 262 END DO 258 263 ! … … 273 278 sshn(:,:) = sshb(:,:) ! set now ssh to the before value 274 279 ! 275 IF( nn_rstssh /= 0 ) THEN 276 nn_rstssh = 0 ! hand-made initilization of ssh 280 IF( nn_rstssh /= 0 ) THEN 281 nn_rstssh = 0 ! hand-made initilization of ssh 277 282 CALL ctl_warn( 'istate_eel: force nn_rstssh = 0' ) 278 283 ENDIF … … 326 331 !!---------------------------------------------------------------------- 327 332 !! *** ROUTINE istate_gyre *** 328 !! 333 !! 329 334 !! ** Purpose : Initialization of the dynamics and tracers for GYRE 330 335 !! configuration (double gyre with rotated domain) … … 351 356 & * (-TANH( (500-fsdept(ji,jj,jk)) / 150 ) + 1) / 2 & 352 357 & + ( 15. * ( 1. - TANH( (fsdept(ji,jj,jk)-50.) / 1500.) ) & 353 & - 1.4 * TANH((fsdept(ji,jj,jk)-100.) / 100.) & 354 & + 7. * (1500. - fsdept(ji,jj,jk)) / 1500. ) & 358 & - 1.4 * TANH((fsdept(ji,jj,jk)-100.) / 100.) & 359 & + 7. * (1500. - fsdept(ji,jj,jk)) / 1500. ) & 355 360 & * (-TANH( (fsdept(ji,jj,jk) - 500) / 150) + 1) / 2 356 361 tsn(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk) … … 363 368 & + 0.2 * TANH( (fsdept(ji,jj,jk) - 35. ) / 100. ) & 364 369 & + 0.2 * TANH( (fsdept(ji,jj,jk) - 1000.) / 5000.) ) & 365 & * (-TANH((fsdept(ji,jj,jk) - 500) / 150) + 1) / 2 370 & * (-TANH((fsdept(ji,jj,jk) - 500) / 150) + 1) / 2 366 371 tsn(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk) 367 372 tsb(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal) … … 379 384 ! ---------------------- 380 385 CALL iom_open ( 'data_tem', inum ) 381 CALL iom_get ( inum, jpdom_data, 'votemper', tsn(:,:,:,jp_tem) ) 386 CALL iom_get ( inum, jpdom_data, 'votemper', tsn(:,:,:,jp_tem) ) 382 387 CALL iom_close( inum ) 383 388 384 tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * tmask(:,:,:) 389 tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * tmask(:,:,:) 385 390 tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) 386 391 … … 388 393 ! ------------------- 389 394 CALL iom_open ( 'data_sal', inum ) 390 CALL iom_get ( inum, jpdom_data, 'vosaline', tsn(:,:,:,jp_sal) ) 395 CALL iom_get ( inum, jpdom_data, 'vosaline', tsn(:,:,:,jp_sal) ) 391 396 CALL iom_close( inum ) 392 397 393 tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:) 398 tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:) 394 399 tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) 395 400 … … 412 417 !! ** Purpose : Compute the geostrophic velocities from (tn,sn) fields 413 418 !! 414 !! ** Method : Using the hydrostatic hypothesis the now hydrostatic 419 !! ** Method : Using the hydrostatic hypothesis the now hydrostatic 415 420 !! pressure is computed by integrating the in-situ density from the 416 421 !! surface to the bottom. … … 429 434 CALL wrk_alloc( jpi, jpj, jpk, zprn) 430 435 ! 431 IF(lwp) WRITE(numout,*) 436 IF(lwp) WRITE(numout,*) 432 437 IF(lwp) WRITE(numout,*) 'istate_uvg : Start from Geostrophy' 433 438 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' … … 437 442 438 443 zalfg = 0.5 * grav * rau0 439 444 440 445 zprn(:,:,1) = zalfg * fse3w(:,:,1) * ( 1 + rhd(:,:,1) ) ! Surface value 441 446 … … 443 448 zprn(:,:,jk) = zprn(:,:,jk-1) & 444 449 & + zalfg * fse3w(:,:,jk) * ( 2. + rhd(:,:,jk) + rhd(:,:,jk-1) ) 445 END DO 450 END DO 446 451 447 452 ! Compute geostrophic balance … … 485 490 CALL lbc_lnk( un, 'U', -1. ) 486 491 CALL lbc_lnk( vn, 'V', -1. ) 487 492 488 493 ub(:,:,:) = un(:,:,:) 489 494 vb(:,:,:) = vn(:,:,:) 490 495 491 496 ! WARNING !!!!! 492 497 ! after initializing u and v, we need to calculate the initial streamfunction bsf. … … 516 521 un(:,:,:) = ub(:,:,:) 517 522 vn(:,:,:) = vb(:,:,:) 518 523 519 524 ! Compute the divergence and curl 520 525 -
branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
r4600 r5168 77 77 #endif 78 78 USE tamtrj ! Output trajectory, needed for TAM 79 USE asminc 79 80 80 81 IMPLICIT NONE … … 110 111 !! Madec, 2008, internal report, IPSL. 111 112 !!---------------------------------------------------------------------- 112 INTEGER :: istp 113 INTEGER :: istp, jk ! time step index 113 114 !!---------------------------------------------------------------------- 114 115 ! … … 155 156 CALL Agrif_Step( stp ) ! AGRIF: time stepping 156 157 #else 158 !-------------------------------------------------------------! 159 ! This trick ensures a minimum consistency between mixing coef. 160 ! and other variables 161 IF( lk_zdftke .AND. lk_asminc .AND. ln_asmdin .AND. ( istp == nit000) ) THEN 162 CALL tke_avn ! recompute avt, avm, 163 ! avmu, avmv and dissl (approximation) 164 DO jk = nit000, nit000 + 2 165 CALL zdf_tke( jk ) ; 166 END DO 167 168 END IF 169 !-------------------------------------------------------------! 170 157 171 CALL stp( istp ) ! standard time stepping 158 172 #endif … … 322 336 CALL dyn_nept_init ! simplified form of Neptune effect 323 337 CALL flush(numout) 324 338 CALL asm_inc_rea_nam ! to fix properly neuler value 339 ! in istate_init 325 340 CALL istate_init ! ocean initial state (Dynamics and tracers) 326 341
Note: See TracChangeset
for help on using the changeset viewer.
