- Timestamp:
- 2015-11-30T20:55:41+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
r5621 r5956 14 14 15 15 !!---------------------------------------------------------------------- 16 !! 'key_asminc' : Switch on the assimilation increment interface17 !!----------------------------------------------------------------------18 16 !! asm_inc_init : Initialize the increment arrays and IAU weights 19 17 !! calc_date : Compute the calendar date YYYYMMDD on a given step … … 28 26 USE domvvl ! domain: variable volume level 29 27 USE oce ! Dynamics and active tracers defined in memory 30 USE ldfdyn _oce ! ocean dynamics: lateral physics28 USE ldfdyn ! lateral diffusion: eddy viscosity coefficients 31 29 USE eosbn2 ! Equation of state - in situ and potential density 32 30 USE zpshde ! Partial step : Horizontal Derivative … … 56 54 LOGICAL, PUBLIC, PARAMETER :: lk_asminc = .FALSE. !: No assimilation increments 57 55 #endif 58 LOGICAL, PUBLIC :: ln_bkgwri = .FALSE.!: No output of the background state fields59 LOGICAL, PUBLIC :: ln_asmiau = .FALSE.!: No applying forcing with an assimilation increment60 LOGICAL, PUBLIC :: ln_asmdin = .FALSE.!: No direct initialization61 LOGICAL, PUBLIC :: ln_trainc = .FALSE.!: No tracer (T and S) assimilation increments62 LOGICAL, PUBLIC :: ln_dyninc = .FALSE.!: No dynamics (u and v) assimilation increments63 LOGICAL, PUBLIC :: ln_sshinc = .FALSE.!: No sea surface height assimilation increment64 LOGICAL, PUBLIC :: ln_seaiceinc 65 LOGICAL, PUBLIC :: ln_salfix = .FALSE.!: Apply minimum salinity check56 LOGICAL, PUBLIC :: ln_bkgwri !: No output of the background state fields 57 LOGICAL, PUBLIC :: ln_asmiau !: No applying forcing with an assimilation increment 58 LOGICAL, PUBLIC :: ln_asmdin !: No direct initialization 59 LOGICAL, PUBLIC :: ln_trainc !: No tracer (T and S) assimilation increments 60 LOGICAL, PUBLIC :: ln_dyninc !: No dynamics (u and v) assimilation increments 61 LOGICAL, PUBLIC :: ln_sshinc !: No sea surface height assimilation increment 62 LOGICAL, PUBLIC :: ln_seaiceinc !: No sea ice concentration increment 63 LOGICAL, PUBLIC :: ln_salfix !: Apply minimum salinity check 66 64 LOGICAL, PUBLIC :: ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing 67 INTEGER, PUBLIC :: nn_divdmp 65 INTEGER, PUBLIC :: nn_divdmp !: Apply divergence damping filter nn_divdmp times 68 66 69 67 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: t_bkg , s_bkg !: Background temperature and salinity … … 90 88 !! * Substitutions 91 89 # include "domzgr_substitute.h90" 92 # include "ldfdyn_substitute.h90"93 90 # include "vectopt_loop_substitute.h90" 94 91 !!---------------------------------------------------------------------- … … 139 136 ! Read Namelist nam_asminc : assimilation increment interface 140 137 !----------------------------------------------------------------------- 141 ln_seaiceinc = .FALSE.138 ln_seaiceinc = .FALSE. 142 139 ln_temnofreeze = .FALSE. 143 140 … … 428 425 429 426 IF ( ln_dyninc .AND. nn_divdmp > 0 ) THEN 430 431 CALL wrk_alloc( jpi,jpj,hdiv)432 433 DO 434 427 ! 428 CALL wrk_alloc( jpi,jpj, hdiv ) 429 ! 430 DO jt = 1, nn_divdmp 431 ! 435 432 DO jk = 1, jpkm1 436 437 433 hdiv(:,:) = 0._wp 438 439 434 DO jj = 2, jpjm1 440 435 DO ji = fs_2, fs_jpim1 ! vector opt. … … 444 439 + e1v(ji ,jj ) * fse3v(ji ,jj ,jk) * v_bkginc(ji ,jj ,jk) & 445 440 - e1v(ji ,jj-1) * fse3v(ji ,jj-1,jk) * v_bkginc(ji ,jj-1,jk) ) & 446 / ( e1 t(ji,jj) *e2t(ji,jj) * fse3t(ji,jj,jk) )441 / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 447 442 END DO 448 443 END DO 449 450 444 CALL lbc_lnk( hdiv, 'T', 1. ) ! lateral boundary cond. (no sign change) 451 445 ! 452 446 DO jj = 2, jpjm1 453 447 DO ji = fs_2, fs_jpim1 ! vector opt. 454 u_bkginc(ji,jj,jk) = u_bkginc(ji,jj,jk) + 0.2_wp * ( e1 t(ji+1,jj)*e2t(ji+1,jj) * hdiv(ji+1,jj) &455 - e1t(ji ,jj)*e2t(ji ,jj) * hdiv(ji ,jj) ) &456 /e1u(ji,jj) * umask(ji,jj,jk)457 v_bkginc(ji,jj,jk) = v_bkginc(ji,jj,jk) + 0.2_wp * ( e1 t(ji,jj+1)*e2t(ji,jj+1) * hdiv(ji,jj+1) &458 - e1t(ji,jj )*e2t(ji,jj ) * hdiv(ji,jj ) ) &459 /e2v(ji,jj) * vmask(ji,jj,jk)448 u_bkginc(ji,jj,jk) = u_bkginc(ji,jj,jk) + 0.2_wp * ( e1e2t(ji+1,jj) * hdiv(ji+1,jj) & 449 & - e1e2t(ji ,jj) * hdiv(ji ,jj) ) & 450 & * r1_e1u(ji,jj) * umask(ji,jj,jk) 451 v_bkginc(ji,jj,jk) = v_bkginc(ji,jj,jk) + 0.2_wp * ( e1e2t(ji,jj+1) * hdiv(ji,jj+1) & 452 & - e1e2t(ji,jj ) * hdiv(ji,jj ) ) & 453 & * r1_e2v(ji,jj) * vmask(ji,jj,jk) 460 454 END DO 461 455 END DO 462 463 456 END DO 464 457 ! 465 458 END DO 466 467 CALL wrk_dealloc( jpi,jpj,hdiv)468 459 ! 460 CALL wrk_dealloc( jpi,jpj, hdiv ) 461 ! 469 462 ENDIF 470 471 472 463 473 464 !----------------------------------------------------------------------- … … 476 467 477 468 IF ( ln_asmdin ) THEN 478 469 ! 479 470 ALLOCATE( t_bkg(jpi,jpj,jpk) ) 480 471 ALLOCATE( s_bkg(jpi,jpj,jpk) ) … … 482 473 ALLOCATE( v_bkg(jpi,jpj,jpk) ) 483 474 ALLOCATE( ssh_bkg(jpi,jpj) ) 484 485 t_bkg(:,:,:) = 0. 0486 s_bkg(:,:,:) = 0. 0487 u_bkg(:,:,:) = 0. 0488 v_bkg(:,:,:) = 0. 0489 ssh_bkg(:,:) = 0. 0490 475 ! 476 t_bkg(:,:,:) = 0._wp 477 s_bkg(:,:,:) = 0._wp 478 u_bkg(:,:,:) = 0._wp 479 v_bkg(:,:,:) = 0._wp 480 ssh_bkg(:,:) = 0._wp 481 ! 491 482 !-------------------------------------------------------------------- 492 483 ! Read from file the background state at analysis time 493 484 !-------------------------------------------------------------------- 494 485 ! 495 486 CALL iom_open( c_asmdin, inum ) 496 487 ! 497 488 CALL iom_get( inum, 'rdastp', zdate_bkg ) 498 489 ! 499 490 IF(lwp) THEN 500 491 WRITE(numout,*) 501 WRITE(numout,*) 'asm_inc_init : Assimilation background state valid at : ', & 502 & NINT( zdate_bkg ) 492 WRITE(numout,*) 'asm_inc_init : Assimilation background state valid at : ', NINT( zdate_bkg ) 503 493 WRITE(numout,*) '~~~~~~~~~~~~' 504 494 ENDIF 505 495 ! 506 496 IF ( NINT( zdate_bkg ) /= iitdin_date ) & 507 497 & CALL ctl_warn( ' Validity time of assimilation background state does', & 508 498 & ' not agree with Direct Initialization time' ) 509 499 ! 510 500 IF ( ln_trainc ) THEN 511 501 CALL iom_get( inum, jpdom_autoglo, 'tn', t_bkg ) … … 514 504 s_bkg(:,:,:) = s_bkg(:,:,:) * tmask(:,:,:) 515 505 ENDIF 516 506 ! 517 507 IF ( ln_dyninc ) THEN 518 508 CALL iom_get( inum, jpdom_autoglo, 'un', u_bkg ) … … 521 511 v_bkg(:,:,:) = v_bkg(:,:,:) * vmask(:,:,:) 522 512 ENDIF 523 513 ! 524 514 IF ( ln_sshinc ) THEN 525 515 CALL iom_get( inum, jpdom_autoglo, 'sshn', ssh_bkg ) 526 516 ssh_bkg(:,:) = ssh_bkg(:,:) * tmask(:,:,1) 527 517 ENDIF 528 518 ! 529 519 CALL iom_close( inum ) 530 520 ! 531 521 ENDIF 532 522 ! … … 574 564 ! If kt = kit000 - 1 then set the date to the restart date 575 565 IF ( kt == kit000 - 1 ) THEN 576 577 566 kdate = ndastp 578 567 RETURN 579 580 568 ENDIF 581 569 … … 646 634 !! ** Action : 647 635 !!---------------------------------------------------------------------- 648 INTEGER, INTENT(IN) :: kt! Current time step649 ! 650 INTEGER :: ji,jj,jk651 INTEGER :: it636 INTEGER, INTENT(IN) :: kt ! Current time step 637 ! 638 INTEGER :: ji, jj, jk 639 INTEGER :: it 652 640 REAL(wp) :: zincwgt ! IAU weight for current time step 653 641 REAL (wp), DIMENSION(jpi,jpj,jpk) :: fzptnz ! 3d freezing point values 654 642 !!---------------------------------------------------------------------- 655 643 ! 656 644 ! freezing point calculation taken from oc_fz_pt (but calculated for all depths) 657 645 ! used to prevent the applied increments taking the temperature below the local freezing point 658 659 646 DO jk = 1, jpkm1 660 647 CALL eos_fzp( tsn(:,:,jk,jp_sal), fzptnz(:,:,jk), fsdept(:,:,jk) ) 661 648 END DO 662 663 IF ( ln_asmiau ) THEN 664 665 !-------------------------------------------------------------------- 666 ! Incremental Analysis Updating 667 !-------------------------------------------------------------------- 668 649 ! 650 ! !-------------------------------------- 651 IF ( ln_asmiau ) THEN ! Incremental Analysis Updating 652 ! !-------------------------------------- 653 ! 669 654 IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN 670 655 ! 671 656 it = kt - nit000 + 1 672 657 zincwgt = wgtiau(it) / rdt ! IAU weight for the current time step 673 658 ! 674 659 IF(lwp) THEN 675 660 WRITE(numout,*) … … 677 662 WRITE(numout,*) '~~~~~~~~~~~~' 678 663 ENDIF 679 664 ! 680 665 ! Update the tracer tendencies 681 666 DO jk = 1, jpkm1 … … 700 685 ENDIF 701 686 END DO 702 703 ENDIF 704 687 ! 688 ENDIF 689 ! 705 690 IF ( kt == nitiaufin_r + 1 ) THEN ! For bias crcn to work 706 691 DEALLOCATE( t_bkginc ) 707 692 DEALLOCATE( s_bkginc ) 708 693 ENDIF 709 710 711 ELSEIF ( ln_asmdin ) THEN 712 713 !-------------------------------------------------------------------- 714 ! Direct Initialization 715 !-------------------------------------------------------------------- 716 694 ! !-------------------------------------- 695 ELSEIF ( ln_asmdin ) THEN ! Direct Initialization 696 ! !-------------------------------------- 697 ! 717 698 IF ( kt == nitdin_r ) THEN 718 699 ! 719 700 neuler = 0 ! Force Euler forward step 720 701 ! 721 702 ! Initialize the now fields with the background + increment 722 703 IF (ln_temnofreeze) THEN … … 745 726 !!gm 746 727 747 748 728 IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav) & 749 729 & CALL zps_hde ( kt, jpts, tsb, gtsu, gtsv, & ! Partial steps: before horizontal gradient … … 766 746 ENDIF 767 747 ! Perhaps the following call should be in step 768 IF 748 IF ( ln_seaiceinc ) CALL seaice_asm_inc ( kt ) ! apply sea ice concentration increment 769 749 ! 770 750 END SUBROUTINE tra_asm_inc … … 787 767 REAL(wp) :: zincwgt ! IAU weight for current time step 788 768 !!---------------------------------------------------------------------- 789 790 IF ( ln_asmiau ) THEN 791 792 !-------------------------------------------------------------------- 793 ! Incremental Analysis Updating 794 !-------------------------------------------------------------------- 795 769 ! 770 ! !-------------------------------------------- 771 IF ( ln_asmiau ) THEN ! Incremental Analysis Updating 772 ! !-------------------------------------------- 773 ! 796 774 IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN 797 775 ! 798 776 it = kt - nit000 + 1 799 777 zincwgt = wgtiau(it) / rdt ! IAU weight for the current time step 800 778 ! 801 779 IF(lwp) THEN 802 780 WRITE(numout,*) 803 WRITE(numout,*) 'dyn_asm_inc : Dynamics IAU at time step = ', & 804 & kt,' with IAU weight = ', wgtiau(it) 781 WRITE(numout,*) 'dyn_asm_inc : Dynamics IAU at time step = ', kt,' with IAU weight = ', wgtiau(it) 805 782 WRITE(numout,*) '~~~~~~~~~~~~' 806 783 ENDIF 807 784 ! 808 785 ! Update the dynamic tendencies 809 786 DO jk = 1, jpkm1 … … 811 788 va(:,:,jk) = va(:,:,jk) + v_bkginc(:,:,jk) * zincwgt 812 789 END DO 813 790 ! 814 791 IF ( kt == nitiaufin_r ) THEN 815 792 DEALLOCATE( u_bkginc ) 816 793 DEALLOCATE( v_bkginc ) 817 794 ENDIF 818 819 ENDIF 820 821 ELSEIF ( ln_asmdin ) THEN 822 823 !-------------------------------------------------------------------- 824 ! Direct Initialization 825 !-------------------------------------------------------------------- 826 795 ! 796 ENDIF 797 ! !----------------------------------------- 798 ELSEIF ( ln_asmdin ) THEN ! Direct Initialization 799 ! !----------------------------------------- 800 ! 827 801 IF ( kt == nitdin_r ) THEN 828 802 ! 829 803 neuler = 0 ! Force Euler forward step 830 804 ! 831 805 ! Initialize the now fields with the background + increment 832 806 un(:,:,:) = u_bkg(:,:,:) + u_bkginc(:,:,:) 833 807 vn(:,:,:) = v_bkg(:,:,:) + v_bkginc(:,:,:) 834 808 ! 835 809 ub(:,:,:) = un(:,:,:) ! Update before fields 836 810 vb(:,:,:) = vn(:,:,:) 837 811 ! 838 812 DEALLOCATE( u_bkg ) 839 813 DEALLOCATE( v_bkg ) … … 863 837 REAL(wp) :: zincwgt ! IAU weight for current time step 864 838 !!---------------------------------------------------------------------- 865 866 IF ( ln_asmiau ) THEN 867 868 !-------------------------------------------------------------------- 869 ! Incremental Analysis Updating 870 !-------------------------------------------------------------------- 871 839 ! 840 ! !----------------------------------------- 841 IF ( ln_asmiau ) THEN ! Incremental Analysis Updating 842 ! !----------------------------------------- 843 ! 872 844 IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN 873 845 ! 874 846 it = kt - nit000 + 1 875 847 zincwgt = wgtiau(it) / rdt ! IAU weight for the current time step 876 848 ! 877 849 IF(lwp) THEN 878 850 WRITE(numout,*) … … 881 853 WRITE(numout,*) '~~~~~~~~~~~~' 882 854 ENDIF 883 855 ! 884 856 ! Save the tendency associated with the IAU weighted SSH increment 885 857 ! (applied in dynspg.*) … … 890 862 DEALLOCATE( ssh_bkginc ) 891 863 ENDIF 892 893 ENDIF 894 895 ELSEIF ( ln_asmdin ) THEN 896 897 !-------------------------------------------------------------------- 898 ! Direct Initialization 899 !-------------------------------------------------------------------- 900 864 ! 865 ENDIF 866 ! !----------------------------------------- 867 ELSEIF ( ln_asmdin ) THEN ! Direct Initialization 868 ! !----------------------------------------- 869 ! 901 870 IF ( kt == nitdin_r ) THEN 902 903 neuler = 0 ! Force Euler forward step 904 905 ! Initialize the now fields the background + increment 906 sshn(:,:) = ssh_bkg(:,:) + ssh_bkginc(:,:) 907 908 ! Update before fields 909 sshb(:,:) = sshn(:,:) 910 871 ! 872 neuler = 0 ! Force Euler forward step 873 ! 874 sshn(:,:) = ssh_bkg(:,:) + ssh_bkginc(:,:) ! Initialize the now fields the background + increment 875 ! 876 sshb(:,:) = sshn(:,:) ! Update before fields 877 ! 911 878 IF( lk_vvl ) THEN 912 879 DO jk = 1, jpk … … 914 881 END DO 915 882 ENDIF 916 883 ! 917 884 DEALLOCATE( ssh_bkg ) 918 885 DEALLOCATE( ssh_bkginc ) 919 886 ! 920 887 ENDIF 921 888 ! … … 936 903 !! 937 904 !!---------------------------------------------------------------------- 938 IMPLICIT NONE 939 ! 940 INTEGER, INTENT(in) :: kt ! Current time step 905 INTEGER, INTENT(in) :: kt ! Current time step 941 906 INTEGER, INTENT(in), OPTIONAL :: kindic ! flag for disabling the deallocation 942 907 ! … … 948 913 #endif 949 914 !!---------------------------------------------------------------------- 950 951 IF ( ln_asmiau ) THEN 952 953 !-------------------------------------------------------------------- 954 ! Incremental Analysis Updating 955 !-------------------------------------------------------------------- 956 915 ! 916 ! !----------------------------------------- 917 IF ( ln_asmiau ) THEN ! Incremental Analysis Updating 918 ! !----------------------------------------- 919 ! 957 920 IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN 958 921 ! 959 922 it = kt - nit000 + 1 960 923 zincwgt = wgtiau(it) ! IAU weight for the current time step 961 924 ! note this is not a tendency so should not be divided by rdt (as with the tracer and other increments) 962 925 ! 963 926 IF(lwp) THEN 964 927 WRITE(numout,*) 965 WRITE(numout,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', & 966 & kt,' with IAU weight = ', wgtiau(it) 928 WRITE(numout,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', kt,' with IAU weight = ', wgtiau(it) 967 929 WRITE(numout,*) '~~~~~~~~~~~~' 968 930 ENDIF 969 931 ! 970 932 ! Sea-ice : LIM-3 case (to add) 971 933 ! 972 934 #if defined key_lim2 973 935 ! Sea-ice : LIM-2 case … … 1007 969 1008 970 #if defined key_cice && defined key_asminc 1009 ! Sea-ice : CICE case. Zero ice increment tendency into CICE 1010 ndaice_da(:,:) = 0.0_wp 1011 #endif 1012 1013 ENDIF 1014 1015 ELSEIF ( ln_asmdin ) THEN 1016 1017 !-------------------------------------------------------------------- 1018 ! Direct Initialization 1019 !-------------------------------------------------------------------- 1020 971 ndaice_da(:,:) = 0._wp ! Sea-ice : CICE case. Zero ice increment tendency into CICE 972 #endif 973 974 ENDIF 975 ! !----------------------------------------- 976 ELSEIF ( ln_asmdin ) THEN ! Direct Initialization 977 ! !----------------------------------------- 978 ! 1021 979 IF ( kt == nitdin_r ) THEN 1022 980 ! 1023 981 neuler = 0 ! Force Euler forward step 1024 982 ! 1025 983 ! Sea-ice : LIM-3 case (to add) 1026 984 ! 1027 985 #if defined key_lim2 1028 986 ! Sea-ice : LIM-2 case. … … 1040 998 zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt 1041 999 ELSEWHERE 1042 zhicifinc(:,:) = 0. 0_wp1000 zhicifinc(:,:) = 0._wp 1043 1001 END WHERE 1044 1002 ! … … 1049 1007 ! seaice salinity balancing (to add) 1050 1008 #endif 1051 1009 ! 1052 1010 #if defined key_cice && defined key_asminc 1053 1011 ! Sea-ice : CICE case. Pass ice increment tendency into CICE 1054 1012 ndaice_da(:,:) = seaice_bkginc(:,:) / rdt 1055 1013 #endif 1056 IF ( .NOT. PRESENT(kindic) ) THEN1057 DEALLOCATE( seaice_bkginc )1058 END IF1059 1014 IF ( .NOT. PRESENT(kindic) ) THEN 1015 DEALLOCATE( seaice_bkginc ) 1016 END IF 1017 ! 1060 1018 ELSE 1061 1019 ! 1062 1020 #if defined key_cice && defined key_asminc 1063 ! Sea-ice : CICE case. Zero ice increment tendency into CICE1064 ndaice_da(:,:) = 0.0_wp 1065 #endif 1066 1021 ndaice_da(:,:) = 0._wp ! Sea-ice : CICE case. Zero ice increment tendency into CICE 1022 1023 #endif 1024 ! 1067 1025 ENDIF 1068 1026 … … 1141 1099 ! 1142 1100 !#endif 1143 1101 ! 1144 1102 ENDIF 1145 1103 ! 1146 1104 END SUBROUTINE seaice_asm_inc 1147 1105
Note: See TracChangeset
for help on using the changeset viewer.
