- Timestamp:
- 2015-10-06T13:40:42+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
r5737 r5777 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. … … 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 448 u_bkginc(ji,jj,jk) = u_bkginc(ji,jj,jk) + 0.2_wp * ( e1e2t(ji+1,jj) * hdiv(ji+1,jj) & 455 456 449 & - e1e2t(ji ,jj) * hdiv(ji ,jj) ) & 450 & * r1_e1u(ji,jj) * umask(ji,jj,jk) 457 451 v_bkginc(ji,jj,jk) = v_bkginc(ji,jj,jk) + 0.2_wp * ( e1e2t(ji,jj+1) * hdiv(ji,jj+1) & 458 459 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 IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav) & 749 & CALL zps_hde ( kt, jpts, tsb, gtsu, gtsv, & ! Partial steps: before horizontal gradient 750 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 751 IF( ln_zps .AND. .NOT. lk_c1d .AND. ln_isfcav) & 752 & CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv, & ! Partial steps for top cell (ISF) 753 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & 754 & gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) ! of t, s, rd at the last ocean level 728 IF( ln_zps .AND. .NOT. lk_c1d ) THEN ! Partial steps: before horizontal gradient 729 IF(ln_isfcav) THEN ! ocean cavities: top and bottom cells (ISF) 730 CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv, gtui, gtvi, & 731 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & 732 & grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) 733 ELSE ! no ocean cavities: bottom cells 734 CALL zps_hde ( kt, jpts, tsb, gtsu, gtsv, & ! 735 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 736 ENDIF 737 ENDIF 755 738 756 739 #if defined key_zdfkpp … … 758 741 !!gm fabien CALL eos( tsn, rhd ) ! Compute rhd 759 742 #endif 760 743 ! 761 744 DEALLOCATE( t_bkginc ) 762 745 DEALLOCATE( s_bkginc ) … … 767 750 ENDIF 768 751 ! Perhaps the following call should be in step 769 IF 752 IF ( ln_seaiceinc ) CALL seaice_asm_inc ( kt ) ! apply sea ice concentration increment 770 753 ! 771 754 END SUBROUTINE tra_asm_inc … … 788 771 REAL(wp) :: zincwgt ! IAU weight for current time step 789 772 !!---------------------------------------------------------------------- 790 791 IF ( ln_asmiau ) THEN 792 793 !-------------------------------------------------------------------- 794 ! Incremental Analysis Updating 795 !-------------------------------------------------------------------- 796 773 ! 774 ! !-------------------------------------------- 775 IF ( ln_asmiau ) THEN ! Incremental Analysis Updating 776 ! !-------------------------------------------- 777 ! 797 778 IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN 798 779 ! 799 780 it = kt - nit000 + 1 800 781 zincwgt = wgtiau(it) / rdt ! IAU weight for the current time step 801 782 ! 802 783 IF(lwp) THEN 803 784 WRITE(numout,*) 804 WRITE(numout,*) 'dyn_asm_inc : Dynamics IAU at time step = ', & 805 & kt,' with IAU weight = ', wgtiau(it) 785 WRITE(numout,*) 'dyn_asm_inc : Dynamics IAU at time step = ', kt,' with IAU weight = ', wgtiau(it) 806 786 WRITE(numout,*) '~~~~~~~~~~~~' 807 787 ENDIF 808 788 ! 809 789 ! Update the dynamic tendencies 810 790 DO jk = 1, jpkm1 … … 812 792 va(:,:,jk) = va(:,:,jk) + v_bkginc(:,:,jk) * zincwgt 813 793 END DO 814 794 ! 815 795 IF ( kt == nitiaufin_r ) THEN 816 796 DEALLOCATE( u_bkginc ) 817 797 DEALLOCATE( v_bkginc ) 818 798 ENDIF 819 820 ENDIF 821 822 ELSEIF ( ln_asmdin ) THEN 823 824 !-------------------------------------------------------------------- 825 ! Direct Initialization 826 !-------------------------------------------------------------------- 827 799 ! 800 ENDIF 801 ! !----------------------------------------- 802 ELSEIF ( ln_asmdin ) THEN ! Direct Initialization 803 ! !----------------------------------------- 804 ! 828 805 IF ( kt == nitdin_r ) THEN 829 806 ! 830 807 neuler = 0 ! Force Euler forward step 831 808 ! 832 809 ! Initialize the now fields with the background + increment 833 810 un(:,:,:) = u_bkg(:,:,:) + u_bkginc(:,:,:) 834 811 vn(:,:,:) = v_bkg(:,:,:) + v_bkginc(:,:,:) 835 812 ! 836 813 ub(:,:,:) = un(:,:,:) ! Update before fields 837 814 vb(:,:,:) = vn(:,:,:) 838 815 ! 839 816 DEALLOCATE( u_bkg ) 840 817 DEALLOCATE( v_bkg ) … … 864 841 REAL(wp) :: zincwgt ! IAU weight for current time step 865 842 !!---------------------------------------------------------------------- 866 867 IF ( ln_asmiau ) THEN 868 869 !-------------------------------------------------------------------- 870 ! Incremental Analysis Updating 871 !-------------------------------------------------------------------- 872 843 ! 844 ! !----------------------------------------- 845 IF ( ln_asmiau ) THEN ! Incremental Analysis Updating 846 ! !----------------------------------------- 847 ! 873 848 IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN 874 849 ! 875 850 it = kt - nit000 + 1 876 851 zincwgt = wgtiau(it) / rdt ! IAU weight for the current time step 877 852 ! 878 853 IF(lwp) THEN 879 854 WRITE(numout,*) … … 882 857 WRITE(numout,*) '~~~~~~~~~~~~' 883 858 ENDIF 884 859 ! 885 860 ! Save the tendency associated with the IAU weighted SSH increment 886 861 ! (applied in dynspg.*) … … 891 866 DEALLOCATE( ssh_bkginc ) 892 867 ENDIF 893 894 ENDIF 895 896 ELSEIF ( ln_asmdin ) THEN 897 898 !-------------------------------------------------------------------- 899 ! Direct Initialization 900 !-------------------------------------------------------------------- 901 868 ! 869 ENDIF 870 ! !----------------------------------------- 871 ELSEIF ( ln_asmdin ) THEN ! Direct Initialization 872 ! !----------------------------------------- 873 ! 902 874 IF ( kt == nitdin_r ) THEN 903 904 neuler = 0 ! Force Euler forward step 905 906 ! Initialize the now fields the background + increment 907 sshn(:,:) = ssh_bkg(:,:) + ssh_bkginc(:,:) 908 909 ! Update before fields 910 sshb(:,:) = sshn(:,:) 911 875 ! 876 neuler = 0 ! Force Euler forward step 877 ! 878 sshn(:,:) = ssh_bkg(:,:) + ssh_bkginc(:,:) ! Initialize the now fields the background + increment 879 ! 880 sshb(:,:) = sshn(:,:) ! Update before fields 881 ! 912 882 IF( lk_vvl ) THEN 913 883 DO jk = 1, jpk … … 915 885 END DO 916 886 ENDIF 917 887 ! 918 888 DEALLOCATE( ssh_bkg ) 919 889 DEALLOCATE( ssh_bkginc ) 920 890 ! 921 891 ENDIF 922 892 ! … … 937 907 !! 938 908 !!---------------------------------------------------------------------- 939 IMPLICIT NONE 940 ! 941 INTEGER, INTENT(in) :: kt ! Current time step 909 INTEGER, INTENT(in) :: kt ! Current time step 942 910 INTEGER, INTENT(in), OPTIONAL :: kindic ! flag for disabling the deallocation 943 911 ! … … 949 917 #endif 950 918 !!---------------------------------------------------------------------- 951 952 IF ( ln_asmiau ) THEN 953 954 !-------------------------------------------------------------------- 955 ! Incremental Analysis Updating 956 !-------------------------------------------------------------------- 957 919 ! 920 ! !----------------------------------------- 921 IF ( ln_asmiau ) THEN ! Incremental Analysis Updating 922 ! !----------------------------------------- 923 ! 958 924 IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN 959 925 ! 960 926 it = kt - nit000 + 1 961 927 zincwgt = wgtiau(it) ! IAU weight for the current time step 962 928 ! note this is not a tendency so should not be divided by rdt (as with the tracer and other increments) 963 929 ! 964 930 IF(lwp) THEN 965 931 WRITE(numout,*) 966 WRITE(numout,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', & 967 & kt,' with IAU weight = ', wgtiau(it) 932 WRITE(numout,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', kt,' with IAU weight = ', wgtiau(it) 968 933 WRITE(numout,*) '~~~~~~~~~~~~' 969 934 ENDIF 970 935 ! 971 936 ! Sea-ice : LIM-3 case (to add) 972 937 ! 973 938 #if defined key_lim2 974 939 ! Sea-ice : LIM-2 case … … 1008 973 1009 974 #if defined key_cice && defined key_asminc 1010 ! Sea-ice : CICE case. Zero ice increment tendency into CICE 1011 ndaice_da(:,:) = 0.0_wp 1012 #endif 1013 1014 ENDIF 1015 1016 ELSEIF ( ln_asmdin ) THEN 1017 1018 !-------------------------------------------------------------------- 1019 ! Direct Initialization 1020 !-------------------------------------------------------------------- 1021 975 ndaice_da(:,:) = 0._wp ! Sea-ice : CICE case. Zero ice increment tendency into CICE 976 #endif 977 978 ENDIF 979 ! !----------------------------------------- 980 ELSEIF ( ln_asmdin ) THEN ! Direct Initialization 981 ! !----------------------------------------- 982 ! 1022 983 IF ( kt == nitdin_r ) THEN 1023 984 ! 1024 985 neuler = 0 ! Force Euler forward step 1025 986 ! 1026 987 ! Sea-ice : LIM-3 case (to add) 1027 988 ! 1028 989 #if defined key_lim2 1029 990 ! Sea-ice : LIM-2 case. … … 1041 1002 zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt 1042 1003 ELSEWHERE 1043 zhicifinc(:,:) = 0. 0_wp1004 zhicifinc(:,:) = 0._wp 1044 1005 END WHERE 1045 1006 ! … … 1050 1011 ! seaice salinity balancing (to add) 1051 1012 #endif 1052 1013 ! 1053 1014 #if defined key_cice && defined key_asminc 1054 1015 ! Sea-ice : CICE case. Pass ice increment tendency into CICE 1055 1016 ndaice_da(:,:) = seaice_bkginc(:,:) / rdt 1056 1017 #endif 1057 IF ( .NOT. PRESENT(kindic) ) THEN1058 DEALLOCATE( seaice_bkginc )1059 END IF1060 1018 IF ( .NOT. PRESENT(kindic) ) THEN 1019 DEALLOCATE( seaice_bkginc ) 1020 END IF 1021 ! 1061 1022 ELSE 1062 1023 ! 1063 1024 #if defined key_cice && defined key_asminc 1064 ! Sea-ice : CICE case. Zero ice increment tendency into CICE1065 ndaice_da(:,:) = 0.0_wp 1066 #endif 1067 1025 ndaice_da(:,:) = 0._wp ! Sea-ice : CICE case. Zero ice increment tendency into CICE 1026 1027 #endif 1028 ! 1068 1029 ENDIF 1069 1030 … … 1142 1103 ! 1143 1104 !#endif 1144 1105 ! 1145 1106 ENDIF 1146 1107 ! 1147 1108 END SUBROUTINE seaice_asm_inc 1148 1109
Note: See TracChangeset
for help on using the changeset viewer.