Changeset 11223
- Timestamp:
- 2019-07-05T20:53:14+02:00 (5 years ago)
- Location:
- NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization
- Files:
-
- 15 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/cfgs/SHARED/namelist_ref
r10808 r11223 609 609 ln_vol = .false. ! total volume correction (see nn_volctl parameter) 610 610 nn_volctl = 1 ! = 0, the total water flux across open boundaries is zero 611 nb_jpk_bdy = -1 ! number of levels in the bdy data (set < 0 if consistent with planned run)612 611 / 613 612 !----------------------------------------------------------------------- -
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/cfgs/SPITZ12/EXPREF/namelist_cfg
r10911 r11223 177 177 nn_rimwidth = 1 ! width of the relaxation zone 178 178 ln_vol = .false. ! total volume correction (see nn_volctl parameter) 179 nb_jpk_bdy = -1 ! number of levels in the bdy data (set < 0 if consistent with planned run)180 179 / 181 180 !----------------------------------------------------------------------- -
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icevar.F90
r10994 r11223 73 73 PUBLIC ice_var_zapneg 74 74 PUBLIC ice_var_roundoff 75 PUBLIC ice_var_itd76 PUBLIC ice_var_itd277 75 PUBLIC ice_var_bv 78 76 PUBLIC ice_var_enthalpy 79 77 PUBLIC ice_var_sshdyn 78 PUBLIC ice_var_itd 79 80 INTERFACE ice_var_itd 81 MODULE PROCEDURE ice_var_itd_1c1c, ice_var_itd_Nc1c, ice_var_itd_1cMc, ice_var_itd_NcMc 82 END INTERFACE 80 83 81 84 !!---------------------------------------------------------------------- … … 656 659 END SUBROUTINE ice_var_roundoff 657 660 661 662 SUBROUTINE ice_var_bv 663 !!------------------------------------------------------------------- 664 !! *** ROUTINE ice_var_bv *** 665 !! 666 !! ** Purpose : computes mean brine volume (%) in sea ice 667 !! 668 !! ** Method : e = - 0.054 * S (ppt) / T (C) 669 !! 670 !! References : Vancoppenolle et al., JGR, 2007 671 !!------------------------------------------------------------------- 672 INTEGER :: ji, jj, jk, jl ! dummy loop indices 673 !!------------------------------------------------------------------- 674 ! 675 !!gm I prefere to use WHERE / ELSEWHERE to set it to zero only where needed <<<=== to be done 676 !! instead of setting everything to zero as just below 677 bv_i (:,:,:) = 0._wp 678 DO jl = 1, jpl 679 DO jk = 1, nlay_i 680 WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 ) 681 bv_i(:,:,jl) = bv_i(:,:,jl) - rTmlt * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 ) 682 END WHERE 683 END DO 684 END DO 685 WHERE( vt_i(:,:) > epsi20 ) ; bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:) 686 ELSEWHERE ; bvm_i(:,:) = 0._wp 687 END WHERE 688 ! 689 END SUBROUTINE ice_var_bv 690 691 692 SUBROUTINE ice_var_enthalpy 693 !!------------------------------------------------------------------- 694 !! *** ROUTINE ice_var_enthalpy *** 695 !! 696 !! ** Purpose : Computes sea ice energy of melting q_i (J.m-3) from temperature 697 !! 698 !! ** Method : Formula (Bitz and Lipscomb, 1999) 699 !!------------------------------------------------------------------- 700 INTEGER :: ji, jk ! dummy loop indices 701 REAL(wp) :: ztmelts ! local scalar 702 !!------------------------------------------------------------------- 703 ! 704 DO jk = 1, nlay_i ! Sea ice energy of melting 705 DO ji = 1, npti 706 ztmelts = - rTmlt * sz_i_1d(ji,jk) 707 t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point => likely conservation issue 708 ! (sometimes zdf scheme produces abnormally high temperatures) 709 e_i_1d(ji,jk) = rhoi * ( rcpi * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) ) & 710 & + rLfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) ) & 711 & - rcp * ztmelts ) 712 END DO 713 END DO 714 DO jk = 1, nlay_s ! Snow energy of melting 715 DO ji = 1, npti 716 e_s_1d(ji,jk) = rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus ) 717 END DO 718 END DO 719 ! 720 END SUBROUTINE ice_var_enthalpy 721 658 722 659 SUBROUTINE ice_var_itd( zhti, zhts, zati, zh_i, zh_s, za_i ) 660 !!------------------------------------------------------------------- 661 !! *** ROUTINE ice_var_itd *** 662 !! 663 !! ** Purpose : converting 1-cat ice to multiple ice categories 723 FUNCTION ice_var_sshdyn(pssh, psnwice_mass, psnwice_mass_b) 724 !!--------------------------------------------------------------------- 725 !! *** ROUTINE ice_var_sshdyn *** 726 !! 727 !! ** Purpose : compute the equivalent ssh in lead when sea ice is embedded 728 !! 729 !! ** Method : ssh_lead = ssh + (Mice + Msnow) / rau0 730 !! 731 !! ** Reference : Jean-Michel Campin, John Marshall, David Ferreira, 732 !! Sea ice-ocean coupling using a rescaled vertical coordinate z*, 733 !! Ocean Modelling, Volume 24, Issues 1-2, 2008 734 !!---------------------------------------------------------------------- 735 ! 736 ! input 737 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pssh !: ssh [m] 738 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass !: mass of snow and ice at current ice time step [Kg/m2] 739 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass_b !: mass of snow and ice at previous ice time step [Kg/m2] 740 ! 741 ! output 742 REAL(wp), DIMENSION(jpi,jpj) :: ice_var_sshdyn ! equivalent ssh in lead [m] 743 ! 744 ! temporary 745 REAL(wp) :: zintn, zintb ! time interpolation weights [] 746 REAL(wp), DIMENSION(jpi,jpj) :: zsnwiceload ! snow and ice load [m] 747 ! 748 ! compute ice load used to define the equivalent ssh in lead 749 IF( ln_ice_embd ) THEN 750 ! 751 ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1} 752 ! = (1/nn_fsbc)^2 * {SUM[n] , n=0,nn_fsbc-1} 753 zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp 754 ! 755 ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1} 756 ! = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1}) 757 zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp 758 ! 759 zsnwiceload(:,:) = ( zintn * psnwice_mass(:,:) + zintb * psnwice_mass_b(:,:) ) * r1_rau0 760 ! 761 ELSE 762 zsnwiceload(:,:) = 0.0_wp 763 ENDIF 764 ! compute equivalent ssh in lead 765 ice_var_sshdyn(:,:) = pssh(:,:) + zsnwiceload(:,:) 766 ! 767 END FUNCTION ice_var_sshdyn 768 769 770 !!------------------------------------------------------------------- 771 !! *** INTERFACE ice_var_itd *** 772 !! 773 !! ** Purpose : converting N-cat ice to jpl ice categories 774 !!------------------------------------------------------------------- 775 SUBROUTINE ice_var_itd_1c1c( zhti, zhts, zati, zh_i, zh_s, za_i ) 776 !!------------------------------------------------------------------- 777 !! ** Purpose : converting 1-cat ice to 1 ice category 778 !!------------------------------------------------------------------- 779 REAL(wp), DIMENSION(:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables 780 REAL(wp), DIMENSION(:), INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables 781 !!------------------------------------------------------------------- 782 zh_i(:) = zhti(:) 783 zh_s(:) = zhts(:) 784 za_i(:) = zati(:) 785 END SUBROUTINE ice_var_itd_1c1c 786 787 SUBROUTINE ice_var_itd_Nc1c( zhti, zhts, zati, zh_i, zh_s, za_i ) 788 !!------------------------------------------------------------------- 789 !! ** Purpose : converting N-cat ice to 1 ice category 790 !!------------------------------------------------------------------- 791 REAL(wp), DIMENSION(:,:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables 792 REAL(wp), DIMENSION(:) , INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables 793 !!------------------------------------------------------------------- 794 ! 795 za_i(:) = SUM( zati(:,:), dim=2 ) 796 ! 797 WHERE( za_i(:) /= 0._wp ) 798 zh_i(:) = SUM( zhti(:,:) * zati(:,:), dim=2 ) / za_i(:) 799 zh_s(:) = SUM( zhts(:,:) * zati(:,:), dim=2 ) / za_i(:) 800 ELSEWHERE 801 zh_i(:) = 0._wp 802 zh_s(:) = 0._wp 803 END WHERE 804 ! 805 END SUBROUTINE ice_var_itd_Nc1c 806 807 SUBROUTINE ice_var_itd_1cMc( zhti, zhts, zati, zh_i, zh_s, za_i ) 808 !!------------------------------------------------------------------- 809 !! 810 !! ** Purpose : converting 1-cat ice to jpl ice categories 664 811 !! 665 812 !! ice thickness distribution follows a gaussian law … … 705 852 zh_s(1:idim,1:jpl) = 0._wp 706 853 za_i(1:idim,1:jpl) = 0._wp 854 ! 855 IF( jpl == 1 ) THEN 856 CALL ice_var_itd_1c1c( zhti, zhts, zati, zh_i(:,1), zh_s(:,1), za_i(:,1) ) 857 RETURN 858 ENDIF 707 859 ! 708 860 DO ji = 1, idim … … 801 953 END DO 802 954 ! 803 END SUBROUTINE ice_var_itd 804 805 806 SUBROUTINE ice_var_itd2( zhti, zhts, zati, zh_i, zh_s, za_i ) 807 !!------------------------------------------------------------------- 808 !! *** ROUTINE ice_var_itd2 *** 955 END SUBROUTINE ice_var_itd_1cMc 956 957 SUBROUTINE ice_var_itd_NcMc( zhti, zhts, zati, zh_i, zh_s, za_i ) 958 !!------------------------------------------------------------------- 809 959 !! 810 960 !! ** Purpose : converting N-cat ice to jpl ice categories … … 845 995 idim = SIZE( zhti, 1 ) 846 996 icat = SIZE( zhti, 2 ) 847 ! 848 ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) ) ! allocate arrays 849 ALLOCATE( jlmin(idim), jlmax(idim) ) 850 851 ! --- initialize output fields to 0 --- ! 852 zh_i(1:idim,1:jpl) = 0._wp 853 zh_s(1:idim,1:jpl) = 0._wp 854 za_i(1:idim,1:jpl) = 0._wp 855 ! 856 ! --- fill the categories --- ! 857 ! find where cat-input = cat-output and fill cat-output fields 858 jlmax(:) = 0 859 jlmin(:) = 999 860 jlfil(:,:) = 0 861 DO jl1 = 1, jpl 862 DO jl2 = 1, icat 997 ! ! ---------------------- ! 998 IF( icat == jpl ) THEN ! input cat = output cat ! 999 ! ! ---------------------- ! 1000 zh_i(:,:) = zhti(:,:) 1001 zh_s(:,:) = zhts(:,:) 1002 za_i(:,:) = zati(:,:) 1003 ! ! ---------------------- ! 1004 ELSEIF( icat == 1 ) THEN ! specific case if N = 1 ! 1005 ! ! ---------------------- ! 1006 ! 1007 CALL ice_var_itd_1cMc( zhti(:,1), zhts(:,1), zati(:,1), zh_i, zh_s, za_i ) 1008 ! 1009 ! ! ---------------------- ! 1010 ELSEIF( jpl == 1 ) THEN ! specific case if M = 1 ! 1011 ! ! ---------------------- ! 1012 ! 1013 CALL ice_var_itd_Nc1c( zhti, zhts, zati, zh_i(:,1), zh_s(:,1), za_i(:,1) ) 1014 ! 1015 ! ! ----------------------- ! 1016 ELSE ! input cat /= output cat ! 1017 ! ! ----------------------- ! 1018 1019 ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) ) ! allocate arrays 1020 ALLOCATE( jlmin(idim), jlmax(idim) ) 1021 1022 ! --- initialize output fields to 0 --- ! 1023 zh_i(1:idim,1:jpl) = 0._wp 1024 zh_s(1:idim,1:jpl) = 0._wp 1025 za_i(1:idim,1:jpl) = 0._wp 1026 ! 1027 ! --- fill the categories --- ! 1028 ! find where cat-input = cat-output and fill cat-output fields 1029 jlmax(:) = 0 1030 jlmin(:) = 999 1031 jlfil(:,:) = 0 1032 DO jl1 = 1, jpl 1033 DO jl2 = 1, icat 1034 DO ji = 1, idim 1035 IF( hi_max(jl1-1) <= zhti(ji,jl2) .AND. hi_max(jl1) > zhti(ji,jl2) ) THEN 1036 ! fill the right category 1037 zh_i(ji,jl1) = zhti(ji,jl2) 1038 zh_s(ji,jl1) = zhts(ji,jl2) 1039 za_i(ji,jl1) = zati(ji,jl2) 1040 ! record categories that are filled 1041 jlmax(ji) = MAX( jlmax(ji), jl1 ) 1042 jlmin(ji) = MIN( jlmin(ji), jl1 ) 1043 jlfil(ji,jl1) = jl1 1044 ENDIF 1045 END DO 1046 END DO 1047 END DO 1048 ! 1049 ! --- fill the gaps between categories --- ! 1050 ! transfer from categories filled at the previous step to the empty ones in between 1051 DO ji = 1, idim 1052 jl1 = jlmin(ji) 1053 jl2 = jlmax(ji) 1054 IF( jl1 > 1 ) THEN 1055 ! fill the lower cat (jl1-1) 1056 za_i(ji,jl1-1) = ztrans * za_i(ji,jl1) 1057 zh_i(ji,jl1-1) = hi_mean(jl1-1) 1058 ! remove from cat jl1 1059 za_i(ji,jl1 ) = ( 1._wp - ztrans ) * za_i(ji,jl1) 1060 ENDIF 1061 IF( jl2 < jpl ) THEN 1062 ! fill the upper cat (jl2+1) 1063 za_i(ji,jl2+1) = ztrans * za_i(ji,jl2) 1064 zh_i(ji,jl2+1) = hi_mean(jl2+1) 1065 ! remove from cat jl2 1066 za_i(ji,jl2 ) = ( 1._wp - ztrans ) * za_i(ji,jl2) 1067 ENDIF 1068 END DO 1069 ! 1070 jlfil2(:,:) = jlfil(:,:) 1071 ! fill categories from low to high 1072 DO jl = 2, jpl-1 863 1073 DO ji = 1, idim 864 IF( hi_max(jl1-1) <= zhti(ji,jl2) .AND. hi_max(jl1) > zhti(ji,jl2) ) THEN 865 ! fill the right category 866 zh_i(ji,jl1) = zhti(ji,jl2) 867 zh_s(ji,jl1) = zhts(ji,jl2) 868 za_i(ji,jl1) = zati(ji,jl2) 869 ! record categories that are filled 870 jlmax(ji) = MAX( jlmax(ji), jl1 ) 871 jlmin(ji) = MIN( jlmin(ji), jl1 ) 872 jlfil(ji,jl1) = jl1 1074 IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN 1075 ! fill high 1076 za_i(ji,jl) = ztrans * za_i(ji,jl-1) 1077 zh_i(ji,jl) = hi_mean(jl) 1078 jlfil(ji,jl) = jl 1079 ! remove low 1080 za_i(ji,jl-1) = ( 1._wp - ztrans ) * za_i(ji,jl-1) 873 1081 ENDIF 874 1082 END DO 875 1083 END DO 876 END DO 877 ! 878 ! --- fill the gaps between categories --- ! 879 ! transfer from categories filled at the previous step to the empty ones in between 880 DO ji = 1, idim 881 jl1 = jlmin(ji) 882 jl2 = jlmax(ji) 883 IF( jl1 > 1 ) THEN 884 ! fill the lower cat (jl1-1) 885 za_i(ji,jl1-1) = ztrans * za_i(ji,jl1) 886 zh_i(ji,jl1-1) = hi_mean(jl1-1) 887 ! remove from cat jl1 888 za_i(ji,jl1 ) = ( 1._wp - ztrans ) * za_i(ji,jl1) 889 ENDIF 890 IF( jl2 < jpl ) THEN 891 ! fill the upper cat (jl2+1) 892 za_i(ji,jl2+1) = ztrans * za_i(ji,jl2) 893 zh_i(ji,jl2+1) = hi_mean(jl2+1) 894 ! remove from cat jl2 895 za_i(ji,jl2 ) = ( 1._wp - ztrans ) * za_i(ji,jl2) 896 ENDIF 897 END DO 898 ! 899 jlfil2(:,:) = jlfil(:,:) 900 ! fill categories from low to high 901 DO jl = 2, jpl-1 902 DO ji = 1, idim 903 IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN 904 ! fill high 905 za_i(ji,jl) = ztrans * za_i(ji,jl-1) 906 zh_i(ji,jl) = hi_mean(jl) 907 jlfil(ji,jl) = jl 908 ! remove low 909 za_i(ji,jl-1) = ( 1._wp - ztrans ) * za_i(ji,jl-1) 910 ENDIF 911 END DO 912 END DO 913 ! 914 ! fill categories from high to low 915 DO jl = jpl-1, 2, -1 916 DO ji = 1, idim 917 IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN 918 ! fill low 919 za_i(ji,jl) = za_i(ji,jl) + ztrans * za_i(ji,jl+1) 920 zh_i(ji,jl) = hi_mean(jl) 921 jlfil2(ji,jl) = jl 922 ! remove high 923 za_i(ji,jl+1) = ( 1._wp - ztrans ) * za_i(ji,jl+1) 924 ENDIF 925 END DO 926 END DO 927 ! 928 DEALLOCATE( jlfil, jlfil2 ) ! deallocate arrays 929 DEALLOCATE( jlmin, jlmax ) 930 ! 931 END SUBROUTINE ice_var_itd2 932 933 934 SUBROUTINE ice_var_bv 935 !!------------------------------------------------------------------- 936 !! *** ROUTINE ice_var_bv *** 937 !! 938 !! ** Purpose : computes mean brine volume (%) in sea ice 939 !! 940 !! ** Method : e = - 0.054 * S (ppt) / T (C) 941 !! 942 !! References : Vancoppenolle et al., JGR, 2007 943 !!------------------------------------------------------------------- 944 INTEGER :: ji, jj, jk, jl ! dummy loop indices 945 !!------------------------------------------------------------------- 946 ! 947 !!gm I prefere to use WHERE / ELSEWHERE to set it to zero only where needed <<<=== to be done 948 !! instead of setting everything to zero as just below 949 bv_i (:,:,:) = 0._wp 950 DO jl = 1, jpl 951 DO jk = 1, nlay_i 952 WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 ) 953 bv_i(:,:,jl) = bv_i(:,:,jl) - rTmlt * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 ) 954 END WHERE 955 END DO 956 END DO 957 WHERE( vt_i(:,:) > epsi20 ) ; bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:) 958 ELSEWHERE ; bvm_i(:,:) = 0._wp 959 END WHERE 960 ! 961 END SUBROUTINE ice_var_bv 962 963 964 SUBROUTINE ice_var_enthalpy 965 !!------------------------------------------------------------------- 966 !! *** ROUTINE ice_var_enthalpy *** 967 !! 968 !! ** Purpose : Computes sea ice energy of melting q_i (J.m-3) from temperature 969 !! 970 !! ** Method : Formula (Bitz and Lipscomb, 1999) 971 !!------------------------------------------------------------------- 972 INTEGER :: ji, jk ! dummy loop indices 973 REAL(wp) :: ztmelts ! local scalar 974 !!------------------------------------------------------------------- 975 ! 976 DO jk = 1, nlay_i ! Sea ice energy of melting 977 DO ji = 1, npti 978 ztmelts = - rTmlt * sz_i_1d(ji,jk) 979 t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point => likely conservation issue 980 ! (sometimes zdf scheme produces abnormally high temperatures) 981 e_i_1d(ji,jk) = rhoi * ( rcpi * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) ) & 982 & + rLfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) ) & 983 & - rcp * ztmelts ) 984 END DO 985 END DO 986 DO jk = 1, nlay_s ! Snow energy of melting 987 DO ji = 1, npti 988 e_s_1d(ji,jk) = rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus ) 989 END DO 990 END DO 991 ! 992 END SUBROUTINE ice_var_enthalpy 993 994 FUNCTION ice_var_sshdyn(pssh, psnwice_mass, psnwice_mass_b) 995 !!--------------------------------------------------------------------- 996 !! *** ROUTINE ice_var_sshdyn *** 997 !! 998 !! ** Purpose : compute the equivalent ssh in lead when sea ice is embedded 999 !! 1000 !! ** Method : ssh_lead = ssh + (Mice + Msnow) / rau0 1001 !! 1002 !! ** Reference : Jean-Michel Campin, John Marshall, David Ferreira, 1003 !! Sea ice-ocean coupling using a rescaled vertical coordinate z*, 1004 !! Ocean Modelling, Volume 24, Issues 1-2, 2008 1005 !!---------------------------------------------------------------------- 1006 ! 1007 ! input 1008 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pssh !: ssh [m] 1009 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass !: mass of snow and ice at current ice time step [Kg/m2] 1010 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass_b !: mass of snow and ice at previous ice time step [Kg/m2] 1011 ! 1012 ! output 1013 REAL(wp), DIMENSION(jpi,jpj) :: ice_var_sshdyn ! equivalent ssh in lead [m] 1014 ! 1015 ! temporary 1016 REAL(wp) :: zintn, zintb ! time interpolation weights [] 1017 REAL(wp), DIMENSION(jpi,jpj) :: zsnwiceload ! snow and ice load [m] 1018 ! 1019 ! compute ice load used to define the equivalent ssh in lead 1020 IF( ln_ice_embd ) THEN 1021 ! 1022 ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1} 1023 ! = (1/nn_fsbc)^2 * {SUM[n] , n=0,nn_fsbc-1} 1024 zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp 1025 ! 1026 ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1} 1027 ! = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1}) 1028 zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp 1029 ! 1030 zsnwiceload(:,:) = ( zintn * psnwice_mass(:,:) + zintb * psnwice_mass_b(:,:) ) * r1_rau0 1031 ! 1032 ELSE 1033 zsnwiceload(:,:) = 0.0_wp 1084 ! 1085 ! fill categories from high to low 1086 DO jl = jpl-1, 2, -1 1087 DO ji = 1, idim 1088 IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN 1089 ! fill low 1090 za_i(ji,jl) = za_i(ji,jl) + ztrans * za_i(ji,jl+1) 1091 zh_i(ji,jl) = hi_mean(jl) 1092 jlfil2(ji,jl) = jl 1093 ! remove high 1094 za_i(ji,jl+1) = ( 1._wp - ztrans ) * za_i(ji,jl+1) 1095 ENDIF 1096 END DO 1097 END DO 1098 ! 1099 DEALLOCATE( jlfil, jlfil2 ) ! deallocate arrays 1100 DEALLOCATE( jlmin, jlmax ) 1101 ! 1034 1102 ENDIF 1035 ! compute equivalent ssh in lead 1036 ice_var_sshdyn(:,:) = pssh(:,:) + zsnwiceload(:,:) 1037 ! 1038 END FUNCTION ice_var_sshdyn 1039 1103 ! 1104 END SUBROUTINE ice_var_itd_NcMc 1040 1105 1041 1106 #else -
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdy_oce.F90
r11191 r11223 42 42 TYPE, PUBLIC :: OBC_DATA !: Storage for external data 43 43 INTEGER , DIMENSION(2) :: nread 44 LOGICAL :: ll_ssh 45 LOGICAL :: ll_u2d 46 LOGICAL :: ll_v2d 47 LOGICAL :: ll_u3d 48 LOGICAL :: ll_v3d 49 LOGICAL :: ll_tem 50 LOGICAL :: ll_sal 51 LOGICAL :: ll_fvl 44 LOGICAL :: lneed_ssh 45 LOGICAL :: lneed_dyn2d 46 LOGICAL :: lneed_dyn3d 47 LOGICAL :: lneed_tra 48 LOGICAL :: lneed_ice 52 49 REAL(wp), POINTER, DIMENSION(:) :: ssh 53 50 REAL(wp), POINTER, DIMENSION(:) :: u2d … … 57 54 REAL(wp), POINTER, DIMENSION(:,:) :: tem 58 55 REAL(wp), POINTER, DIMENSION(:,:) :: sal 59 #if defined key_si3 60 LOGICAL :: ll_a_i 61 LOGICAL :: ll_h_i 62 LOGICAL :: ll_h_s 63 REAL(wp), POINTER, DIMENSION(:,:) :: a_i !: now ice leads fraction climatology 64 REAL(wp), POINTER, DIMENSION(:,:) :: h_i !: Now ice thickness climatology 65 REAL(wp), POINTER, DIMENSION(:,:) :: h_s !: now snow thickness 66 #endif 56 REAL(wp), POINTER, DIMENSION(:,:) :: a_i !: now ice leads fraction climatology 57 REAL(wp), POINTER, DIMENSION(:,:) :: h_i !: Now ice thickness climatology 58 REAL(wp), POINTER, DIMENSION(:,:) :: h_s !: now snow thickness 67 59 #if defined key_top 68 60 CHARACTER(LEN=20) :: cn_obc !: type of boundary condition to apply … … 87 79 ! 88 80 INTEGER :: nb_bdy !: number of open boundary sets 89 INTEGER, DIMENSION(jp_bdy) :: nb_jpk_bdy !: number of levels in the bdy data (set < 0 if consistent with planned run)90 81 INTEGER, DIMENSION(jp_bdy) :: nn_rimwidth !: boundary rim width for Flow Relaxation Scheme 91 82 INTEGER :: nn_volctl !: = 0 the total volume will have the variability of the surface Flux E-P … … 130 121 INTEGER, DIMENSION(jp_bdy) :: nn_dta !: =0 => *all* data is set to initial conditions 131 122 !: =1 => some data to be read in from data files 132 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: dta_global !: workspace for reading in global data arrays (unstr. bdy)133 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: dta_global_z !: workspace for reading in global depth arrays (unstr. bdy)134 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: dta_global_dz !: workspace for reading in global depth arrays (unstr. bdy)135 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: dta_global2 !: workspace for reading in global data arrays (struct. bdy)136 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: dta_global2_z !: workspace for reading in global depth arrays (struct. bdy)137 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: dta_global2_dz !: workspace for reading in global depth arrays (struct. bdy)138 123 !$AGRIF_DO_NOT_TREAT 139 124 TYPE(OBC_INDEX), DIMENSION(jp_bdy), TARGET :: idx_bdy !: bdy indices (local process) -
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdydta.F90
r10951 r11223 43 43 PUBLIC bdy_dta_init ! routine called by nemogcm.F90 44 44 45 INTEGER, ALLOCATABLE, DIMENSION(:) :: nb_bdy_fld ! Number of fields to update for each boundary set. 46 INTEGER :: nb_bdy_fld_sum ! Total number of fields to update for all boundary sets. 47 LOGICAL, DIMENSION(jp_bdy) :: ln_full_vel_array ! =T => full velocities in 3D boundary conditions 48 ! =F => baroclinic velocities in 3D boundary conditions 45 INTEGER , PARAMETER :: jpbdyfld = 10 ! maximum number of files to read 46 INTEGER , PARAMETER :: jp_bdyssh = 1 ! 47 INTEGER , PARAMETER :: jp_bdyu2d = 2 ! 48 INTEGER , PARAMETER :: jp_bdyv2d = 3 ! 49 INTEGER , PARAMETER :: jp_bdyu3d = 4 ! 50 INTEGER , PARAMETER :: jp_bdyv3d = 5 ! 51 INTEGER , PARAMETER :: jp_bdytem = 6 ! 52 INTEGER , PARAMETER :: jp_bdysal = 7 ! 53 INTEGER , PARAMETER :: jp_bdya_i = 8 ! 54 INTEGER , PARAMETER :: jp_bdyh_i = 9 ! 55 INTEGER , PARAMETER :: jp_bdyh_S = 10 ! 56 ! =F => baroclinic velocities in 3D boundary conditions 49 57 !$AGRIF_DO_NOT_TREAT 50 TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(: ), TARGET :: bf! structure of input fields (file informations, fields read)58 TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:,:), TARGET :: bf ! structure of input fields (file informations, fields read) 51 59 !$AGRIF_END_DO_NOT_TREAT 52 TYPE(MAP_POINTER), ALLOCATABLE, DIMENSION(:) :: nbmap_ptr ! array of pointers to nbmap53 54 #if defined key_si355 INTEGER :: nice_cat ! number of categories in the input file56 INTEGER :: jfld_hti, jfld_hts, jfld_ai ! indices of ice thickness, snow thickness and concentration in bf structure57 INTEGER, DIMENSION(jp_bdy) :: jfld_htit, jfld_htst, jfld_ait58 #endif59 60 60 61 !!---------------------------------------------------------------------- … … 65 66 CONTAINS 66 67 67 SUBROUTINE bdy_dta( kt, jit, time_offset )68 SUBROUTINE bdy_dta( kt, kit, kt_offset ) 68 69 !!---------------------------------------------------------------------- 69 70 !! *** SUBROUTINE bdy_dta *** … … 75 76 !!---------------------------------------------------------------------- 76 77 INTEGER, INTENT(in) :: kt ! ocean time-step index 77 INTEGER, INTENT(in), OPTIONAL :: jit ! subcycle time-step index (for timesplitting option)78 INTEGER, INTENT(in), OPTIONAL :: time_offset ! time offset in units of timesteps. NB. if jit78 INTEGER, INTENT(in), OPTIONAL :: kit ! subcycle time-step index (for timesplitting option) 79 INTEGER, INTENT(in), OPTIONAL :: kt_offset ! time offset in units of timesteps. NB. if kit 79 80 ! ! is present then units = subcycle timesteps. 80 ! ! time_offset = 0 => get data at "now" time level81 ! ! time_offset = -1 => get data at "before" time level82 ! ! time_offset = +1 => get data at "after" time level81 ! ! kt_offset = 0 => get data at "now" time level 82 ! ! kt_offset = -1 => get data at "before" time level 83 ! ! kt_offset = +1 => get data at "after" time level 83 84 ! ! etc. 84 85 ! 85 INTEGER :: jbdy, jfld, jstart, jend, ib, jl ! dummy loop indices 86 INTEGER :: ii, ij, ik, igrd ! local integers 87 INTEGER, DIMENSION(jpbgrd) :: ilen1 88 INTEGER, POINTER, DIMENSION(:) :: nblen, nblenrim ! short cuts 89 TYPE(OBC_DATA), POINTER :: dta ! short cut 86 INTEGER :: jbdy, jfld, jstart, jend, ib, jl ! dummy loop indices 87 INTEGER :: ii, ij, ik, igrd, ipl ! local integers 88 INTEGER, DIMENSION(jpbgrd) :: ilen1 89 INTEGER, DIMENSION(:), POINTER :: nblen, nblenrim ! short cuts 90 TYPE(OBC_DATA) , POINTER :: dta_alias ! short cut 91 TYPE(FLD), DIMENSION(:), POINTER :: bf_alias 90 92 !!--------------------------------------------------------------------------- 91 93 ! … … 94 96 ! Initialise data arrays once for all from initial conditions where required 95 97 !--------------------------------------------------------------------------- 96 IF( kt == nit000 .AND. .NOT.PRESENT( jit) ) THEN98 IF( kt == nit000 .AND. .NOT.PRESENT(kit) ) THEN 97 99 98 100 ! Calculate depth-mean currents 99 101 !----------------------------- 100 102 101 103 DO jbdy = 1, nb_bdy 102 104 ! 103 105 nblen => idx_bdy(jbdy)%nblen 104 106 nblenrim => idx_bdy(jbdy)%nblenrim 105 dta => dta_bdy(jbdy)106 107 ! 107 108 IF( nn_dyn2d_dta(jbdy) == 0 ) THEN 108 109 ilen1(:) = nblen(:) 109 IF( dta %ll_ssh ) THEN110 IF( dta_bdy(jbdy)%lneed_ssh ) THEN 110 111 igrd = 1 111 112 DO ib = 1, ilen1(igrd) … … 113 114 ij = idx_bdy(jbdy)%nbj(ib,igrd) 114 115 dta_bdy(jbdy)%ssh(ib) = sshn(ii,ij) * tmask(ii,ij,1) 115 END DO 116 ENDIF 117 IF( dta %ll_u2d) THEN116 END DO 117 ENDIF 118 IF( dta_bdy(jbdy)%lneed_dyn2d) THEN 118 119 igrd = 2 119 120 DO ib = 1, ilen1(igrd) … … 121 122 ij = idx_bdy(jbdy)%nbj(ib,igrd) 122 123 dta_bdy(jbdy)%u2d(ib) = un_b(ii,ij) * umask(ii,ij,1) 123 END DO 124 ENDIF 125 IF( dta%ll_v2d ) THEN 124 END DO 126 125 igrd = 3 127 126 DO ib = 1, ilen1(igrd) … … 129 128 ij = idx_bdy(jbdy)%nbj(ib,igrd) 130 129 dta_bdy(jbdy)%v2d(ib) = vn_b(ii,ij) * vmask(ii,ij,1) 131 END DO 130 END DO 132 131 ENDIF 133 132 ENDIF … … 135 134 IF( nn_dyn3d_dta(jbdy) == 0 ) THEN 136 135 ilen1(:) = nblen(:) 137 IF( dta %ll_u3d ) THEN136 IF( dta_bdy(jbdy)%lneed_dyn3d ) THEN 138 137 igrd = 2 139 138 DO ib = 1, ilen1(igrd) … … 143 142 dta_bdy(jbdy)%u3d(ib,ik) = ( un(ii,ij,ik) - un_b(ii,ij) ) * umask(ii,ij,ik) 144 143 END DO 145 END DO 146 ENDIF 147 IF( dta%ll_v3d ) THEN 144 END DO 148 145 igrd = 3 149 146 DO ib = 1, ilen1(igrd) … … 152 149 ij = idx_bdy(jbdy)%nbj(ib,igrd) 153 150 dta_bdy(jbdy)%v3d(ib,ik) = ( vn(ii,ij,ik) - vn_b(ii,ij) ) * vmask(ii,ij,ik) 154 155 END DO 151 END DO 152 END DO 156 153 ENDIF 157 154 ENDIF … … 159 156 IF( nn_tra_dta(jbdy) == 0 ) THEN 160 157 ilen1(:) = nblen(:) 161 IF( dta %ll_tem) THEN158 IF( dta_bdy(jbdy)%lneed_tra ) THEN 162 159 igrd = 1 163 160 DO ib = 1, ilen1(igrd) … … 165 162 ii = idx_bdy(jbdy)%nbi(ib,igrd) 166 163 ij = idx_bdy(jbdy)%nbj(ib,igrd) 167 dta_bdy(jbdy)%tem(ib,ik) = tsn(ii,ij,ik,jp_tem) * tmask(ii,ij,ik) 164 dta_bdy(jbdy)%tem(ib,ik) = tsn(ii,ij,ik,jp_bdytem) * tmask(ii,ij,ik) 165 dta_bdy(jbdy)%sal(ib,ik) = tsn(ii,ij,ik,jp_bdysal) * tmask(ii,ij,ik) 168 166 END DO 169 END DO 170 ENDIF 171 IF( dta%ll_sal ) THEN 172 igrd = 1 173 DO ib = 1, ilen1(igrd) 174 DO ik = 1, jpkm1 175 ii = idx_bdy(jbdy)%nbi(ib,igrd) 176 ij = idx_bdy(jbdy)%nbj(ib,igrd) 177 dta_bdy(jbdy)%sal(ib,ik) = tsn(ii,ij,ik,jp_sal) * tmask(ii,ij,ik) 178 END DO 179 END DO 167 END DO 180 168 ENDIF 181 169 ENDIF … … 184 172 IF( nn_ice_dta(jbdy) == 0 ) THEN ! set ice to initial values 185 173 ilen1(:) = nblen(:) 186 IF( dta %ll_a_i) THEN174 IF( dta_bdy(jbdy)%lneed_ice ) THEN 187 175 igrd = 1 188 176 DO jl = 1, jpl … … 191 179 ij = idx_bdy(jbdy)%nbj(ib,igrd) 192 180 dta_bdy(jbdy)%a_i (ib,jl) = a_i(ii,ij,jl) * tmask(ii,ij,1) 193 END DO194 END DO195 ENDIF196 IF( dta%ll_h_i ) THEN197 igrd = 1198 DO jl = 1, jpl199 DO ib = 1, ilen1(igrd)200 ii = idx_bdy(jbdy)%nbi(ib,igrd)201 ij = idx_bdy(jbdy)%nbj(ib,igrd)202 181 dta_bdy(jbdy)%h_i (ib,jl) = h_i(ii,ij,jl) * tmask(ii,ij,1) 203 END DO204 END DO205 ENDIF206 IF( dta%ll_h_s ) THEN207 igrd = 1208 DO jl = 1, jpl209 DO ib = 1, ilen1(igrd)210 ii = idx_bdy(jbdy)%nbi(ib,igrd)211 ij = idx_bdy(jbdy)%nbj(ib,igrd)212 182 dta_bdy(jbdy)%h_s (ib,jl) = h_s(ii,ij,jl) * tmask(ii,ij,1) 213 183 END DO … … 222 192 ! update external data from files 223 193 !-------------------------------- 224 225 jstart = 1 226 DO jbdy = 1, nb_bdy 227 dta => dta_bdy(jbdy) 228 IF( nn_dta(jbdy) == 1 ) THEN ! skip this bit if no external data required 229 230 IF( PRESENT(jit) ) THEN 231 ! Update barotropic boundary conditions only 232 ! jit is optional argument for fld_read and bdytide_update 233 IF( cn_dyn2d(jbdy) /= 'none' ) THEN 234 IF( nn_dyn2d_dta(jbdy) == 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays 235 IF( dta%ll_ssh ) dta%ssh(:) = 0._wp 236 IF( dta%ll_u2d ) dta%u2d(:) = 0._wp 237 IF( dta%ll_u3d ) dta%v2d(:) = 0._wp 238 ENDIF 239 IF (cn_tra(jbdy) /= 'runoff') THEN 240 IF( nn_dyn2d_dta(jbdy) == 1 .OR. nn_dyn2d_dta(jbdy) == 3 ) THEN 241 242 jend = jstart + dta%nread(2) - 1 243 IF( ln_full_vel_array(jbdy) ) THEN 244 CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), & 245 & kit=jit, kt_offset=time_offset , jpk_bdy=nb_jpk_bdy(jbdy), & 246 & fvl=ln_full_vel_array(jbdy) ) 247 ELSE 248 CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), & 249 & kit=jit, kt_offset=time_offset ) 250 ENDIF 251 252 ! If full velocities in boundary data then extract barotropic velocities from 3D fields 253 IF( ln_full_vel_array(jbdy) .AND. & 254 & ( nn_dyn2d_dta(jbdy) == 1 .OR. nn_dyn2d_dta(jbdy) == 3 .OR. & 255 & nn_dyn3d_dta(jbdy) == 1 ) )THEN 256 257 igrd = 2 ! zonal velocity 258 dta%u2d(:) = 0._wp 259 DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 260 ii = idx_bdy(jbdy)%nbi(ib,igrd) 261 ij = idx_bdy(jbdy)%nbj(ib,igrd) 262 DO ik = 1, jpkm1 263 dta%u2d(ib) = dta%u2d(ib) & 264 & + e3u_n(ii,ij,ik) * umask(ii,ij,ik) * dta%u3d(ib,ik) 265 END DO 266 dta%u2d(ib) = dta%u2d(ib) * r1_hu_n(ii,ij) 267 END DO 268 igrd = 3 ! meridional velocity 269 dta%v2d(:) = 0._wp 270 DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 271 ii = idx_bdy(jbdy)%nbi(ib,igrd) 272 ij = idx_bdy(jbdy)%nbj(ib,igrd) 273 DO ik = 1, jpkm1 274 dta%v2d(ib) = dta%v2d(ib) & 275 & + e3v_n(ii,ij,ik) * vmask(ii,ij,ik) * dta%v3d(ib,ik) 276 END DO 277 dta%v2d(ib) = dta%v2d(ib) * r1_hv_n(ii,ij) 278 END DO 279 ENDIF 280 ENDIF 281 IF( nn_dyn2d_dta(jbdy) .ge. 2 ) THEN ! update tidal harmonic forcing 282 CALL bdytide_update( kt=kt, idx=idx_bdy(jbdy), dta=dta, td=tides(jbdy), & 283 & jit=jit, time_offset=time_offset ) 284 ENDIF 285 ENDIF 286 ENDIF 287 ELSE 288 IF (cn_tra(jbdy) == 'runoff') then ! runoff condition 289 jend = nb_bdy_fld(jbdy) 290 CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), & 291 & map=nbmap_ptr(jstart:jend), kt_offset=time_offset ) 292 ! 293 igrd = 2 ! zonal velocity 294 DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 295 ii = idx_bdy(jbdy)%nbi(ib,igrd) 296 ij = idx_bdy(jbdy)%nbj(ib,igrd) 297 dta%u2d(ib) = dta%u2d(ib) / ( e2u(ii,ij) * hu_0(ii,ij) ) 298 END DO 299 ! 300 igrd = 3 ! meridional velocity 301 DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 302 ii = idx_bdy(jbdy)%nbi(ib,igrd) 303 ij = idx_bdy(jbdy)%nbj(ib,igrd) 304 dta%v2d(ib) = dta%v2d(ib) / ( e1v(ii,ij) * hv_0(ii,ij) ) 305 END DO 306 ELSE 307 IF( nn_dyn2d_dta(jbdy) == 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays 308 IF( dta%ll_ssh ) dta%ssh(:) = 0._wp 309 IF( dta%ll_u2d ) dta%u2d(:) = 0._wp 310 IF( dta%ll_v2d ) dta%v2d(:) = 0._wp 311 ENDIF 312 IF( dta%nread(1) .gt. 0 ) THEN ! update external data 313 jend = jstart + dta%nread(1) - 1 314 CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), & 315 & map=nbmap_ptr(jstart:jend), kt_offset=time_offset, jpk_bdy=nb_jpk_bdy(jbdy), & 316 & fvl=ln_full_vel_array(jbdy) ) 317 ENDIF 318 ! If full velocities in boundary data then split into barotropic and baroclinic data 319 IF( ln_full_vel_array(jbdy) .and. & 320 & ( nn_dyn2d_dta(jbdy) == 1 .OR. nn_dyn2d_dta(jbdy) == 3 .OR. & 321 & nn_dyn3d_dta(jbdy) == 1 ) ) THEN 322 igrd = 2 ! zonal velocity 323 dta%u2d(:) = 0._wp 324 DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 325 ii = idx_bdy(jbdy)%nbi(ib,igrd) 326 ij = idx_bdy(jbdy)%nbj(ib,igrd) 327 DO ik = 1, jpkm1 328 dta%u2d(ib) = dta%u2d(ib) & 329 & + e3u_n(ii,ij,ik) * umask(ii,ij,ik) * dta%u3d(ib,ik) 330 END DO 331 dta%u2d(ib) = dta%u2d(ib) * r1_hu_n(ii,ij) 332 DO ik = 1, jpkm1 333 dta%u3d(ib,ik) = dta%u3d(ib,ik) - dta%u2d(ib) 334 END DO 335 END DO 336 igrd = 3 ! meridional velocity 337 dta%v2d(:) = 0._wp 338 DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 339 ii = idx_bdy(jbdy)%nbi(ib,igrd) 340 ij = idx_bdy(jbdy)%nbj(ib,igrd) 341 DO ik = 1, jpkm1 342 dta%v2d(ib) = dta%v2d(ib) & 343 & + e3v_n(ii,ij,ik) * vmask(ii,ij,ik) * dta%v3d(ib,ik) 344 END DO 345 dta%v2d(ib) = dta%v2d(ib) * r1_hv_n(ii,ij) 346 DO ik = 1, jpkm1 347 dta%v3d(ib,ik) = dta%v3d(ib,ik) - dta%v2d(ib) 348 END DO 349 END DO 350 ENDIF 351 352 ENDIF 194 195 DO jbdy = 1, nb_bdy 196 197 dta_alias => dta_bdy(jbdy) 198 bf_alias => bf(:,jbdy) 199 200 ! read/update all bdy data 201 ! ------------------------ 202 CALL fld_read( kt, 1, bf_alias, kit = kit, kt_offset = kt_offset ) 203 204 ! apply some corrections in some specific cases... 205 ! -------------------------------------------------- 206 ! 207 ! if runoff condition: change river flow we read (in m3/s) into barotropic velocity (m/s) 208 IF( cn_tra(jbdy) == 'runoff' .AND. TRIM(bf_alias(jp_bdyu2d)%clrootname) /= 'NOT USED' ) THEN ! runoff and we read u/v2d 209 ! 210 igrd = 2 ! zonal flow (m3/s) to barotropic zonal velocity (m/s) 211 DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 212 ii = idx_bdy(jbdy)%nbi(ib,igrd) 213 ij = idx_bdy(jbdy)%nbj(ib,igrd) 214 dta_alias%u2d(ib) = dta_alias%u2d(ib) / ( e2u(ii,ij) * hu_0(ii,ij) ) 215 END DO 216 igrd = 3 ! meridional flow (m3/s) to barotropic meridional velocity (m/s) 217 DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 218 ii = idx_bdy(jbdy)%nbi(ib,igrd) 219 ij = idx_bdy(jbdy)%nbj(ib,igrd) 220 dta_alias%v2d(ib) = dta_alias%v2d(ib) / ( e1v(ii,ij) * hv_0(ii,ij) ) 221 END DO 222 ENDIF 223 224 ! tidal harmonic forcing ONLY: initialise arrays 225 IF( nn_dyn2d_dta(jbdy) == 2 ) THEN ! we did not read ssh, u/v2d 226 IF( dta_alias%lneed_ssh ) dta_alias%ssh(:) = 0._wp 227 IF( dta_alias%lneed_dyn2d ) dta_alias%u2d(:) = 0._wp 228 IF( dta_alias%lneed_dyn2d ) dta_alias%v2d(:) = 0._wp 229 ENDIF 230 231 ! If full velocities in boundary data, then split it into barotropic and baroclinic component 232 IF( bf_alias(jp_bdyu3d)%ltotvel ) THEN ! if we read 3D total velocity (can be true only if u3d was read) 233 ! 234 igrd = 2 ! zonal velocity 235 dta_alias%u2d(:) = 0._wp ! compute barotrope zonal velocity and put it in u2d 236 DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 237 ii = idx_bdy(jbdy)%nbi(ib,igrd) 238 ij = idx_bdy(jbdy)%nbj(ib,igrd) 239 DO ik = 1, jpkm1 240 dta_alias%u2d(ib) = dta_alias%u2d(ib) + e3u_n(ii,ij,ik) * umask(ii,ij,ik) * dta_alias%u3d(ib,ik) 241 END DO 242 dta_alias%u2d(ib) = dta_alias%u2d(ib) * r1_hu_n(ii,ij) 243 DO ik = 1, jpkm1 ! compute barocline zonal velocity and put it in u3d 244 dta_alias%u3d(ib,ik) = dta_alias%u3d(ib,ik) - dta_alias%u2d(ib) 245 END DO 246 END DO 247 igrd = 3 ! meridional velocity 248 dta_alias%v2d(:) = 0._wp ! compute barotrope meridional velocity and put it in v2d 249 DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 250 ii = idx_bdy(jbdy)%nbi(ib,igrd) 251 ij = idx_bdy(jbdy)%nbj(ib,igrd) 252 DO ik = 1, jpkm1 253 dta_alias%v2d(ib) = dta_alias%v2d(ib) + e3v_n(ii,ij,ik) * vmask(ii,ij,ik) * dta_alias%v3d(ib,ik) 254 END DO 255 dta_alias%v2d(ib) = dta_alias%v2d(ib) * r1_hv_n(ii,ij) 256 DO ik = 1, jpkm1 ! compute barocline meridional velocity and put it in v3d 257 dta_alias%v3d(ib,ik) = dta_alias%v3d(ib,ik) - dta_alias%v2d(ib) 258 END DO 259 END DO 260 ENDIF ! ltotvel 261 262 ! update tidal harmonic forcing 263 IF( PRESENT(kit) .AND. nn_dyn2d_dta(jbdy) .GE. 2 ) THEN 264 CALL bdytide_update( kt = kt, idx = idx_bdy(jbdy), dta = dta_alias, td = tides(jbdy), & 265 & kit = kit, kt_offset = kt_offset ) 266 ENDIF 267 268 ! atm surface pressure : add inverted barometer effect to ssh if it was read 269 IF ( ln_apr_obc .AND. TRIM(bf_alias(jp_bdyssh)%clrootname) /= 'NOT USED' ) THEN 270 igrd = 1 271 DO ib = 1, idx_bdy(jbdy)%nblenrim(igrd) ! ssh is used only on the rim 272 ii = idx_bdy(jbdy)%nbi(ib,igrd) 273 ij = idx_bdy(jbdy)%nbj(ib,igrd) 274 dta_alias%ssh(ib) = dta_alias%ssh(ib) + ssh_ib(ii,ij) 275 END DO 276 ENDIF 277 353 278 #if defined key_si3 354 ! convert N-cat fields (input) into jpl-cat (output) 355 IF( cn_ice(jbdy) /= 'none' .AND. nn_ice_dta(jbdy) == 1 ) THEN 356 jfld_hti = jfld_htit(jbdy) 357 jfld_hts = jfld_htst(jbdy) 358 jfld_ai = jfld_ait(jbdy) 359 IF ( jpl /= 1 .AND. nice_cat == 1 ) THEN ! case input cat = 1 360 CALL ice_var_itd ( bf(jfld_hti)%fnow(:,1,1), bf(jfld_hts)%fnow(:,1,1), bf(jfld_ai)%fnow(:,1,1), & 361 & dta_bdy(jbdy)%h_i , dta_bdy(jbdy)%h_s , dta_bdy(jbdy)%a_i ) 362 ELSEIF( jpl /= 1 .AND. nice_cat /= 1 .AND. nice_cat /= jpl ) THEN ! case input cat /=1 and /=jpl 363 CALL ice_var_itd2( bf(jfld_hti)%fnow(:,1,:), bf(jfld_hts)%fnow(:,1,:), bf(jfld_ai)%fnow(:,1,:), & 364 & dta_bdy(jbdy)%h_i , dta_bdy(jbdy)%h_s , dta_bdy(jbdy)%a_i ) 365 ENDIF 366 ENDIF 279 ! ice: convert N-cat fields (input) into jpl-cat (output) 280 IF( dta_alias%lneed_ice ) THEN 281 ipl = SIZE(bf_alias(jp_bdya_i)%fnow, 3) 282 IF( ipl /= jpl ) THEN ! ice: convert N-cat fields (input) into jpl-cat (output) 283 CALL ice_var_itd(bf_alias(jp_bdyh_i)%fnow(:,1,:), bf_alias(jp_bdyh_s)%fnow(:,1,:), bf_alias(jp_bdya_i)%fnow(:,1,:), & 284 & dta_alias%h_i , dta_alias%h_s , dta_alias%a_i ) 285 ENDIF 286 ENDIF 367 287 #endif 368 ENDIF369 jstart = jstart + dta%nread(1)370 ENDIF ! nn_dta(jbdy) = 1371 288 END DO ! jbdy 372 373 IF ( ln_apr_obc ) THEN374 DO jbdy = 1, nb_bdy375 IF (cn_tra(jbdy) /= 'runoff')THEN376 igrd = 1 ! meridional velocity377 DO ib = 1, idx_bdy(jbdy)%nblenrim(igrd)378 ii = idx_bdy(jbdy)%nbi(ib,igrd)379 ij = idx_bdy(jbdy)%nbj(ib,igrd)380 dta_bdy(jbdy)%ssh(ib) = dta_bdy(jbdy)%ssh(ib) + ssh_ib(ii,ij)381 END DO382 ENDIF383 END DO384 ENDIF385 289 386 290 IF ( ln_tide ) THEN 387 291 IF (ln_dynspg_ts) THEN ! Fill temporary arrays with slow-varying bdy data 388 DO jbdy = 1, nb_bdy ! Tidal component added in ts loop389 IF ( nn_dyn2d_dta(jbdy) . ge. 2 ) THEN292 DO jbdy = 1, nb_bdy ! Tidal component added in ts loop 293 IF ( nn_dyn2d_dta(jbdy) .GE. 2 ) THEN 390 294 nblen => idx_bdy(jbdy)%nblen 391 295 nblenrim => idx_bdy(jbdy)%nblenrim 392 IF( cn_dyn2d(jbdy) == 'frs' ) THEN; ilen1(:)=nblen(:) ; ELSE ; ilen1(:)=nblenrim(:) ; ENDIF 393 IF ( dta_bdy(jbdy)%ll_ssh ) dta_bdy_s(jbdy)%ssh(1:ilen1(1)) = dta_bdy(jbdy)%ssh(1:ilen1(1)) 394 IF ( dta_bdy(jbdy)%ll_u2d ) dta_bdy_s(jbdy)%u2d(1:ilen1(2)) = dta_bdy(jbdy)%u2d(1:ilen1(2)) 395 IF ( dta_bdy(jbdy)%ll_v2d ) dta_bdy_s(jbdy)%v2d(1:ilen1(3)) = dta_bdy(jbdy)%v2d(1:ilen1(3)) 396 ENDIF 397 END DO 398 ELSE ! Add tides if not split-explicit free surface else this is done in ts loop 399 ! 400 CALL bdy_dta_tides( kt=kt, time_offset=time_offset ) 401 ENDIF 402 ENDIF 403 404 ! 405 IF( ln_timing ) CALL timing_stop('bdy_dta') 406 ! 407 END SUBROUTINE bdy_dta 296 IF( cn_dyn2d(jbdy) == 'frs' ) THEN ; ilen1(:)=nblen(:) ; ELSE ; ilen1(:)=nblenrim(:) ; ENDIF 297 IF ( dta_bdy(jbdy)%lneed_ssh ) dta_bdy_s(jbdy)%ssh(1:ilen1(1)) = dta_bdy(jbdy)%ssh(1:ilen1(1)) 298 IF ( dta_bdy(jbdy)%lneed_dyn2d ) dta_bdy_s(jbdy)%u2d(1:ilen1(2)) = dta_bdy(jbdy)%u2d(1:ilen1(2)) 299 IF ( dta_bdy(jbdy)%lneed_dyn2d ) dta_bdy_s(jbdy)%v2d(1:ilen1(3)) = dta_bdy(jbdy)%v2d(1:ilen1(3)) 300 ENDIF 301 END DO 302 ELSE ! Add tides if not split-explicit free surface else this is done in ts loop 303 ! 304 CALL bdy_dta_tides( kt=kt, kt_offset=kt_offset ) 305 ENDIF 306 ENDIF 307 ! 308 IF( ln_timing ) CALL timing_stop('bdy_dta') 309 ! 310 END SUBROUTINE bdy_dta 408 311 409 312 … … 418 321 !! 419 322 !!---------------------------------------------------------------------- 420 INTEGER :: jbdy, jfld, jstart, jend, ierror, ios ! Local integers 421 ! 323 INTEGER :: jbdy, jfld ! Local integers 324 INTEGER :: ierror, ios ! 325 ! 326 CHARACTER(len=3) :: cl3 ! 422 327 CHARACTER(len=100) :: cn_dir ! Root directory for location of data files 423 CHARACTER(len=100), DIMENSION(nb_bdy) :: cn_dir_array ! Root directory for location of data files424 CHARACTER(len = 256):: clname ! temporary file name425 328 LOGICAL :: ln_full_vel ! =T => full velocities in 3D boundary data 426 329 ! ! =F => baroclinic velocities in 3D boundary data 427 INTEGER :: ilen_global ! Max length required for global bdy dta arrays 428 INTEGER, ALLOCATABLE, DIMENSION(:) :: ilen1, ilen3 ! size of 1st and 3rd dimensions of local arrays 429 INTEGER, ALLOCATABLE, DIMENSION(:) :: ibdy ! bdy set for a particular jfld 430 INTEGER, ALLOCATABLE, DIMENSION(:) :: igrid ! index for grid type (1,2,3 = T,U,V) 431 INTEGER, POINTER, DIMENSION(:) :: nblen, nblenrim ! short cuts 432 TYPE(OBC_DATA), POINTER :: dta ! short cut 433 #if defined key_si3 434 INTEGER :: kndims ! number of dimensions in an array (i.e. 3 = wo ice cat; 4 = w ice cat) 435 INTEGER, DIMENSION(4) :: kdimsz ! size of dimensions 436 INTEGER :: inum,id1 ! local integer 437 #endif 438 TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) :: blf_i ! array of namelist information structures 439 TYPE(FLD_N) :: bn_tem, bn_sal, bn_u3d, bn_v3d ! 440 TYPE(FLD_N) :: bn_ssh, bn_u2d, bn_v2d ! informations about the fields to be read 441 #if defined key_si3 442 TYPE(FLD_N) :: bn_a_i, bn_h_i, bn_h_s 443 #endif 330 INTEGER :: ipk,ipl ! 331 INTEGER :: idvar ! variable ID 332 INTEGER :: indims ! number of dimensions of the variable 333 INTEGER :: iszdim ! number of dimensions of the variable 334 INTEGER, DIMENSION(4) :: i4dimsz ! size of variable dimensions 335 INTEGER :: igrd ! index for grid type (1,2,3 = T,U,V) 336 LOGICAL :: lluld ! is the variable using the unlimited dimension 337 LOGICAL :: llneed ! 338 LOGICAL :: llread ! 339 TYPE(FLD_N), DIMENSION(1), TARGET :: bn_tem, bn_sal, bn_u3d, bn_v3d ! must be an array to be used with fld_fill 340 TYPE(FLD_N), DIMENSION(1), TARGET :: bn_ssh, bn_u2d, bn_v2d ! informations about the fields to be read 341 TYPE(FLD_N), DIMENSION(1), TARGET :: bn_a_i, bn_h_i, bn_h_s 342 TYPE(FLD_N), DIMENSION(:), POINTER :: bn_alias ! must be an array to be used with fld_fill 343 TYPE(FLD ), DIMENSION(:), POINTER :: bf_alias 344 ! 444 345 NAMELIST/nambdy_dta/ cn_dir, bn_tem, bn_sal, bn_u3d, bn_v3d, bn_ssh, bn_u2d, bn_v2d 445 #if defined key_si3446 346 NAMELIST/nambdy_dta/ bn_a_i, bn_h_i, bn_h_s 447 #endif448 347 NAMELIST/nambdy_dta/ ln_full_vel 449 348 !!--------------------------------------------------------------------------- … … 454 353 IF(lwp) WRITE(numout,*) '' 455 354 456 ! Set nn_dta 457 DO jbdy = 1, nb_bdy 458 nn_dta(jbdy) = MAX( nn_dyn2d_dta (jbdy) & 459 & , nn_dyn3d_dta (jbdy) & 460 & , nn_tra_dta (jbdy) & 461 #if defined key_si3 462 & , nn_ice_dta (jbdy) & 463 #endif 464 ) 465 IF(nn_dta(jbdy) > 1) nn_dta(jbdy) = 1 466 END DO 467 468 ! Work out upper bound of how many fields there are to read in and allocate arrays 469 ! --------------------------------------------------------------------------- 470 ALLOCATE( nb_bdy_fld(nb_bdy) ) 471 nb_bdy_fld(:) = 0 472 DO jbdy = 1, nb_bdy 473 IF( cn_dyn2d(jbdy) /= 'none' .AND. ( nn_dyn2d_dta(jbdy) == 1 .or. nn_dyn2d_dta(jbdy) == 3 ) ) THEN 474 nb_bdy_fld(jbdy) = nb_bdy_fld(jbdy) + 3 475 ENDIF 476 IF( cn_dyn3d(jbdy) /= 'none' .AND. nn_dyn3d_dta(jbdy) == 1 ) THEN 477 nb_bdy_fld(jbdy) = nb_bdy_fld(jbdy) + 2 478 ENDIF 479 IF( cn_tra(jbdy) /= 'none' .AND. nn_tra_dta(jbdy) == 1 ) THEN 480 nb_bdy_fld(jbdy) = nb_bdy_fld(jbdy) + 2 481 ENDIF 482 #if defined key_si3 483 IF( cn_ice(jbdy) /= 'none' .AND. nn_ice_dta(jbdy) == 1 ) THEN 484 nb_bdy_fld(jbdy) = nb_bdy_fld(jbdy) + 3 485 ENDIF 486 #endif 487 IF(lwp) WRITE(numout,*) 'Maximum number of files to open =', nb_bdy_fld(jbdy) 488 END DO 489 490 nb_bdy_fld_sum = SUM( nb_bdy_fld ) 491 492 ALLOCATE( bf(nb_bdy_fld_sum), STAT=ierror ) 355 ALLOCATE( bf(jpbdyfld,nb_bdy), STAT=ierror ) 493 356 IF( ierror > 0 ) THEN 494 357 CALL ctl_stop( 'bdy_dta: unable to allocate bf structure' ) ; RETURN 495 358 ENDIF 496 ALLOCATE( blf_i(nb_bdy_fld_sum), STAT=ierror ) 497 IF( ierror > 0 ) THEN 498 CALL ctl_stop( 'bdy_dta: unable to allocate blf_i structure' ) ; RETURN 499 ENDIF 500 ALLOCATE( nbmap_ptr(nb_bdy_fld_sum), STAT=ierror ) 501 IF( ierror > 0 ) THEN 502 CALL ctl_stop( 'bdy_dta: unable to allocate nbmap_ptr structure' ) ; RETURN 503 ENDIF 504 ALLOCATE( ilen1(nb_bdy_fld_sum), ilen3(nb_bdy_fld_sum) ) 505 ALLOCATE( ibdy(nb_bdy_fld_sum) ) 506 ALLOCATE( igrid(nb_bdy_fld_sum) ) 507 359 bf(:,:)%clrootname = 'NOT USED' ! default definition used as a flag in fld_read to do nothing. 360 bf(:,:)%ltotvel = .FALSE. ! default definition 361 508 362 ! Read namelists 509 363 ! -------------- 510 364 REWIND(numnam_cfg) 511 jfld = 0 512 DO jbdy = 1, nb_bdy 513 IF( nn_dta(jbdy) == 1 ) THEN 514 REWIND(numnam_ref) 515 READ ( numnam_ref, nambdy_dta, IOSTAT = ios, ERR = 901) 516 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy_dta in reference namelist', lwp ) 365 DO jbdy = 1, nb_bdy 366 367 WRITE(ctmp1, '(a,i2)') 'BDY number ', jbdy 368 WRITE(ctmp2, '(a,i2)') 'block nambdy_dta number ', jbdy 369 370 ! There is only one nambdy_dta block in namelist_ref -> use it for each bdy so we do a rewind 371 REWIND(numnam_ref) 372 READ ( numnam_ref, nambdy_dta, IOSTAT = ios, ERR = 901) 373 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy_dta in reference namelist', lwp ) 374 375 ! by-pass nambdy_dta reading if no input data used in this bdy 376 IF( ( dta_bdy(jbdy)%lneed_dyn2d .AND. MOD(nn_dyn2d_dta(jbdy),2) == 1 ) & 377 & .OR. ( dta_bdy(jbdy)%lneed_dyn3d .AND. nn_dyn3d_dta(jbdy) == 1 ) & 378 & .OR. ( dta_bdy(jbdy)%lneed_tra .AND. nn_tra_dta(jbdy) == 1 ) & 379 & .OR. ( dta_bdy(jbdy)%lneed_ice .AND. nn_ice_dta(jbdy) == 1 ) ) THEN 380 ! WARNING: we don't do a rewind here, each bdy reads its own nambdy_dta block one after another 517 381 READ ( numnam_cfg, nambdy_dta, IOSTAT = ios, ERR = 902 ) 518 382 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nambdy_dta in configuration namelist', lwp ) 519 IF(lwm) WRITE( numond, nambdy_dta ) 520 521 cn_dir_array(jbdy) = cn_dir 522 ln_full_vel_array(jbdy) = ln_full_vel 523 524 nblen => idx_bdy(jbdy)%nblen 525 nblenrim => idx_bdy(jbdy)%nblenrim 526 dta => dta_bdy(jbdy) 527 dta%nread(2) = 0 528 529 ! Only read in necessary fields for this set. 530 ! Important that barotropic variables come first. 531 IF( nn_dyn2d_dta(jbdy) == 1 .or. nn_dyn2d_dta(jbdy) == 3 ) THEN 532 533 IF( dta%ll_ssh ) THEN 534 if(lwp) write(numout,*) '++++++ reading in ssh field' 535 jfld = jfld + 1 536 blf_i(jfld) = bn_ssh 537 ibdy(jfld) = jbdy 538 igrid(jfld) = 1 539 ilen1(jfld) = nblen(igrid(jfld)) 540 ilen3(jfld) = 1 541 dta%nread(2) = dta%nread(2) + 1 542 ENDIF 543 544 IF( dta%ll_u2d .and. .not. ln_full_vel_array(jbdy) ) THEN 545 if(lwp) write(numout,*) '++++++ reading in u2d field' 546 jfld = jfld + 1 547 blf_i(jfld) = bn_u2d 548 ibdy(jfld) = jbdy 549 igrid(jfld) = 2 550 ilen1(jfld) = nblen(igrid(jfld)) 551 ilen3(jfld) = 1 552 dta%nread(2) = dta%nread(2) + 1 553 ENDIF 554 555 IF( dta%ll_v2d .and. .not. ln_full_vel_array(jbdy) ) THEN 556 if(lwp) write(numout,*) '++++++ reading in v2d field' 557 jfld = jfld + 1 558 blf_i(jfld) = bn_v2d 559 ibdy(jfld) = jbdy 560 igrid(jfld) = 3 561 ilen1(jfld) = nblen(igrid(jfld)) 562 ilen3(jfld) = 1 563 dta%nread(2) = dta%nread(2) + 1 564 ENDIF 565 566 ENDIF 567 568 ! read 3D velocities if baroclinic velocities require OR if 569 ! barotropic velocities required and ln_full_vel set to .true. 570 IF( nn_dyn3d_dta(jbdy) == 1 .OR. & 571 & ( ln_full_vel_array(jbdy) .AND. ( nn_dyn2d_dta(jbdy) == 1 .OR. nn_dyn2d_dta(jbdy) == 3 ) ) ) THEN 572 573 IF( dta%ll_u3d .OR. ( ln_full_vel_array(jbdy) .and. dta%ll_u2d ) ) THEN 574 if(lwp) write(numout,*) '++++++ reading in u3d field' 575 jfld = jfld + 1 576 blf_i(jfld) = bn_u3d 577 ibdy(jfld) = jbdy 578 igrid(jfld) = 2 579 ilen1(jfld) = nblen(igrid(jfld)) 580 ilen3(jfld) = jpk 581 IF( ln_full_vel_array(jbdy) .and. dta%ll_u2d ) dta%nread(2) = dta%nread(2) + 1 582 ENDIF 583 584 IF( dta%ll_v3d .OR. ( ln_full_vel_array(jbdy) .and. dta%ll_v2d ) ) THEN 585 if(lwp) write(numout,*) '++++++ reading in v3d field' 586 jfld = jfld + 1 587 blf_i(jfld) = bn_v3d 588 ibdy(jfld) = jbdy 589 igrid(jfld) = 3 590 ilen1(jfld) = nblen(igrid(jfld)) 591 ilen3(jfld) = jpk 592 IF( ln_full_vel_array(jbdy) .and. dta%ll_v2d ) dta%nread(2) = dta%nread(2) + 1 593 ENDIF 594 595 ENDIF 596 597 ! temperature and salinity 598 IF( nn_tra_dta(jbdy) == 1 ) THEN 599 600 IF( dta%ll_tem ) THEN 601 if(lwp) write(numout,*) '++++++ reading in tem field' 602 jfld = jfld + 1 603 blf_i(jfld) = bn_tem 604 ibdy(jfld) = jbdy 605 igrid(jfld) = 1 606 ilen1(jfld) = nblen(igrid(jfld)) 607 ilen3(jfld) = jpk 608 ENDIF 609 610 IF( dta%ll_sal ) THEN 611 if(lwp) write(numout,*) '++++++ reading in sal field' 612 jfld = jfld + 1 613 blf_i(jfld) = bn_sal 614 ibdy(jfld) = jbdy 615 igrid(jfld) = 1 616 ilen1(jfld) = nblen(igrid(jfld)) 617 ilen3(jfld) = jpk 618 ENDIF 619 620 ENDIF 621 622 #if defined key_si3 623 ! sea ice 624 IF( nn_ice_dta(jbdy) == 1 ) THEN 625 ! Test for types of ice input (1cat or Xcat) 626 ! Build file name to find dimensions 627 clname=TRIM( cn_dir )//TRIM(bn_a_i%clname) 628 IF( .NOT. bn_a_i%ln_clim ) THEN 629 WRITE(clname, '(a,"_y",i4.4)' ) TRIM( clname ), nyear ! add year 630 IF( bn_a_i%cltype /= 'yearly' ) WRITE(clname, '(a,"m" ,i2.2)' ) TRIM( clname ), nmonth ! add month 631 ELSE 632 IF( bn_a_i%cltype /= 'yearly' ) WRITE(clname, '(a,"_m",i2.2)' ) TRIM( clname ), nmonth ! add month 633 ENDIF 634 IF( bn_a_i%cltype == 'daily' .OR. bn_a_i%cltype(1:4) == 'week' ) & 635 & WRITE(clname, '(a,"d" ,i2.2)' ) TRIM( clname ), nday ! add day 383 IF(lwm) WRITE( numond, nambdy_dta ) 384 ENDIF 385 386 ! get the number of ice categories in bdy data file (use a_i information to do this) 387 ipl = jpl ! default definition 388 IF( dta_bdy(jbdy)%lneed_ice ) THEN ! if we need ice bdy data 389 IF( nn_ice_dta(jbdy) == 1 ) THEN ! if we get ice bdy data from netcdf file 390 CALL fld_fill( bf(jp_bdya_i,jbdy:jbdy), bn_a_i, cn_dir, 'bdy_dta', 'a_i'//' '//ctmp1, ctmp2 ) ! use namelist info 391 CALL fld_clopn( bf(jp_bdya_i,jbdy), nyear, nmonth, nday ) ! not a problem when we call it again after 392 idvar = iom_varid( bf(jp_bdya_i,jbdy)%num, bf(jp_bdya_i,jbdy)%clvar, kndims=indims, kdimsz=i4dimsz, lduld=lluld ) 393 IF( indims == 4 .OR. ( indims == 3 .AND. .NOT. lluld ) ) THEN ; ipl = i4dimsz(3) ! xylt or xyl 394 ELSE ; ipl = 1 ! xy or xyt 395 ENDIF 396 ENDIF 397 ENDIF 398 399 DO jfld = 1, jpbdyfld 400 401 ! ===================== 402 ! ssh 403 ! ===================== 404 IF( jfld == jp_bdyssh ) THEN 405 cl3 = 'ssh' 406 igrd = 1 ! T point 407 ipk = 1 ! surface data 408 llneed = dta_bdy(jbdy)%lneed_ssh ! dta_bdy(jbdy)%ssh will be needed 409 llread = MOD(nn_dyn2d_dta(jbdy),2) == 1 ! get data from NetCDF file 410 bf_alias => bf(jp_bdyssh,jbdy:jbdy) ! alias for ssh structure of bdy number jbdy 411 bn_alias => bn_ssh ! alias for ssh structure of nambdy_dta 412 ENDIF 413 ! ===================== 414 ! dyn2d 415 ! ===================== 416 IF( jfld == jp_bdyu2d ) THEN 417 cl3 = 'u2d' 418 igrd = 2 ! U point 419 ipk = 1 ! surface data 420 llneed = dta_bdy(jbdy)%lneed_dyn2d ! dta_bdy(jbdy)%ssh will be needed 421 llread = .NOT. ln_full_vel .AND. MOD(nn_dyn2d_dta(jbdy),2) == 1 ! don't get u2d from u3d and read NetCDF file 422 bf_alias => bf(jp_bdyu2d,jbdy:jbdy) ! alias for u2d structure of bdy number jbdy 423 bn_alias => bn_u2d ! alias for u2d structure of nambdy_dta 424 ENDIF 425 IF( jfld == jp_bdyv2d ) THEN 426 cl3 = 'v2d' 427 igrd = 3 ! V point 428 ipk = 1 ! surface data 429 llneed = dta_bdy(jbdy)%lneed_dyn2d ! dta_bdy(jbdy)%ssh will be needed 430 llread = .NOT. ln_full_vel .AND. MOD(nn_dyn2d_dta(jbdy),2) == 1 ! don't get v2d from v3d and read NetCDF file 431 bf_alias => bf(jp_bdyv2d,jbdy:jbdy) ! alias for v2d structure of bdy number jbdy 432 bn_alias => bn_v2d ! alias for v2d structure of nambdy_dta 433 ENDIF 434 ! ===================== 435 ! dyn3d 436 ! ===================== 437 IF( jfld == jp_bdyu3d ) THEN 438 cl3 = 'u3d' 439 igrd = 2 ! U point 440 ipk = jpk ! 3d data 441 llneed = dta_bdy(jbdy)%lneed_dyn3d .OR. & ! dta_bdy(jbdy)%u3d will be needed 442 & ( dta_bdy(jbdy)%lneed_dyn2d .AND. ln_full_vel ) ! u3d needed to compute u2d 443 llread = nn_dyn3d_dta(jbdy) == 1 ! get data from NetCDF file 444 bf_alias => bf(jp_bdyu3d,jbdy:jbdy) ! alias for u3d structure of bdy number jbdy 445 bn_alias => bn_u3d ! alias for u3d structure of nambdy_dta 446 ENDIF 447 IF( jfld == jp_bdyv3d ) THEN 448 cl3 = 'v3d' 449 igrd = 3 ! V point 450 ipk = jpk ! 3d data 451 llneed = dta_bdy(jbdy)%lneed_dyn3d .OR. & ! dta_bdy(jbdy)%v3d will be needed 452 & ( dta_bdy(jbdy)%lneed_dyn2d .AND. ln_full_vel ) ! v3d needed to compute v2d 453 llread = nn_dyn3d_dta(jbdy) == 1 ! get data from NetCDF file 454 bf_alias => bf(jp_bdyv3d,jbdy:jbdy) ! alias for v3d structure of bdy number jbdy 455 bn_alias => bn_v3d ! alias for v3d structure of nambdy_dta 456 ENDIF 457 458 ! ===================== 459 ! tra 460 ! ===================== 461 IF( jfld == jp_bdytem ) THEN 462 cl3 = 'tem' 463 igrd = 1 ! T point 464 ipk = jpk ! 3d data 465 llneed = dta_bdy(jbdy)%lneed_tra ! dta_bdy(jbdy)%tem will be needed 466 llread = nn_tra_dta(jbdy) == 1 ! get data from NetCDF file 467 bf_alias => bf(jp_bdytem,jbdy:jbdy) ! alias for ssh structure of bdy number jbdy 468 bn_alias => bn_tem ! alias for ssh structure of nambdy_dta 469 ENDIF 470 IF( jfld == jp_bdysal ) THEN 471 cl3 = 'sal' 472 igrd = 1 ! T point 473 ipk = jpk ! 3d data 474 llneed = dta_bdy(jbdy)%lneed_tra ! dta_bdy(jbdy)%sal will be needed 475 llread = nn_tra_dta(jbdy) == 1 ! get data from NetCDF file 476 bf_alias => bf(jp_bdysal,jbdy:jbdy) ! alias for ssh structure of bdy number jbdy 477 bn_alias => bn_sal ! alias for ssh structure of nambdy_dta 478 ENDIF 479 480 ! ===================== 481 ! ice 482 ! ===================== 483 IF( jfld == jp_bdya_i ) THEN 484 cl3 = 'a_i' 485 igrd = 1 ! T point 486 ipk = ipl ! 487 llneed = dta_bdy(jbdy)%lneed_ice ! dta_bdy(jbdy)%a_i will be needed 488 llread = nn_ice_dta(jbdy) == 1 ! get data from NetCDF file 489 bf_alias => bf(jp_bdya_i,jbdy:jbdy) ! alias for ssh structure of bdy number jbdy 490 bn_alias => bn_a_i ! alias for ssh structure of nambdy_dta 491 ENDIF 492 IF( jfld == jp_bdyh_i ) THEN 493 cl3 = 'h_i' 494 igrd = 1 ! T point 495 ipk = ipl ! 496 llneed = dta_bdy(jbdy)%lneed_ice ! dta_bdy(jbdy)%h_i will be needed 497 llread = nn_ice_dta(jbdy) == 1 ! get data from NetCDF file 498 bf_alias => bf(jp_bdyh_i,jbdy:jbdy) ! alias for ssh structure of bdy number jbdy 499 bn_alias => bn_h_i ! alias for ssh structure of nambdy_dta 500 ENDIF 501 IF( jfld == jp_bdyh_s ) THEN 502 cl3 = 'h_s' 503 igrd = 1 ! T point 504 ipk = ipl ! 505 llneed = dta_bdy(jbdy)%lneed_ice ! dta_bdy(jbdy)%h_s will be needed 506 llread = nn_ice_dta(jbdy) == 1 ! get data from NetCDF file 507 bf_alias => bf(jp_bdyh_s,jbdy:jbdy) ! alias for ssh structure of bdy number jbdy 508 bn_alias => bn_h_s ! alias for ssh structure of nambdy_dta 509 ENDIF 510 511 512 IF( llneed ) THEN ! dta_bdy(jbdy)%xxx will be needed 513 ! ! -> must be associated with an allocated target 514 iszdim = idx_bdy(jbdy)%nblen(igrd) ! length of this bdy on this MPI processus 515 ALLOCATE( bf_alias(1)%fnow( iszdim, 1, ipk ) ) ! allocate the target 636 516 ! 637 CALL iom_open ( clname, inum ) 638 id1 = iom_varid( inum, bn_a_i%clvar, kdimsz=kdimsz, kndims=kndims, ldstop = .FALSE. ) 639 CALL iom_close ( inum ) 640 641 IF ( kndims == 4 ) THEN 642 nice_cat = kdimsz(4) ! Xcat input 643 ELSE 644 nice_cat = 1 ! 1cat input 645 ENDIF 646 ! End test 647 648 IF( dta%ll_a_i ) THEN 649 jfld = jfld + 1 650 blf_i(jfld) = bn_a_i 651 ibdy(jfld) = jbdy 652 igrid(jfld) = 1 653 ilen1(jfld) = nblen(igrid(jfld)) 654 ilen3(jfld) = nice_cat 655 ENDIF 656 657 IF( dta%ll_h_i ) THEN 658 jfld = jfld + 1 659 blf_i(jfld) = bn_h_i 660 ibdy(jfld) = jbdy 661 igrid(jfld) = 1 662 ilen1(jfld) = nblen(igrid(jfld)) 663 ilen3(jfld) = nice_cat 664 ENDIF 665 666 IF( dta%ll_h_s ) THEN 667 jfld = jfld + 1 668 blf_i(jfld) = bn_h_s 669 ibdy(jfld) = jbdy 670 igrid(jfld) = 1 671 ilen1(jfld) = nblen(igrid(jfld)) 672 ilen3(jfld) = nice_cat 673 ENDIF 674 675 ENDIF 676 #endif 677 ! Recalculate field counts 678 !------------------------- 679 IF( jbdy == 1 ) THEN 680 nb_bdy_fld_sum = 0 681 nb_bdy_fld(jbdy) = jfld 682 nb_bdy_fld_sum = jfld 683 ELSE 684 nb_bdy_fld(jbdy) = jfld - nb_bdy_fld_sum 685 nb_bdy_fld_sum = nb_bdy_fld_sum + nb_bdy_fld(jbdy) 686 ENDIF 687 688 dta%nread(1) = nb_bdy_fld(jbdy) 689 690 ENDIF ! nn_dta == 1 691 ENDDO ! jbdy 692 693 DO jfld = 1, nb_bdy_fld_sum 694 ALLOCATE( bf(jfld)%fnow(ilen1(jfld),1,ilen3(jfld)) ) 695 IF( blf_i(jfld)%ln_tint ) ALLOCATE( bf(jfld)%fdta(ilen1(jfld),1,ilen3(jfld),2) ) 696 nbmap_ptr(jfld)%ptr => idx_bdy(ibdy(jfld))%nbmap(:,igrid(jfld)) 697 nbmap_ptr(jfld)%ll_unstruc = ln_coords_file(ibdy(jfld)) 698 ENDDO 699 700 ! fill bf with blf_i and control print 701 !------------------------------------- 702 jstart = 1 703 DO jbdy = 1, nb_bdy 704 jend = jstart - 1 + nb_bdy_fld(jbdy) 705 CALL fld_fill( bf(jstart:jend), blf_i(jstart:jend), cn_dir_array(jbdy), 'bdy_dta', & 706 & 'open boundary conditions', 'nambdy_dta' ) 707 jstart = jend + 1 708 ENDDO 709 710 DO jfld = 1, nb_bdy_fld_sum 711 bf(jfld)%igrd = igrid(jfld) 712 bf(jfld)%ibdy = ibdy(jfld) 713 ENDDO 714 715 ! Initialise local boundary data arrays 716 ! nn_xxx_dta=0 : allocate space - will be filled from initial conditions later 717 ! nn_xxx_dta=1 : point to "fnow" arrays 718 !------------------------------------- 719 720 jfld = 0 721 DO jbdy=1, nb_bdy 722 723 nblen => idx_bdy(jbdy)%nblen 724 dta => dta_bdy(jbdy) 725 726 if(lwp) then 727 write(numout,*) '++++++ dta%ll_ssh = ',dta%ll_ssh 728 write(numout,*) '++++++ dta%ll_u2d = ',dta%ll_u2d 729 write(numout,*) '++++++ dta%ll_v2d = ',dta%ll_v2d 730 write(numout,*) '++++++ dta%ll_u3d = ',dta%ll_u3d 731 write(numout,*) '++++++ dta%ll_v3d = ',dta%ll_v3d 732 write(numout,*) '++++++ dta%ll_tem = ',dta%ll_tem 733 write(numout,*) '++++++ dta%ll_sal = ',dta%ll_sal 734 endif 735 736 IF ( nn_dyn2d_dta(jbdy) == 0 .or. nn_dyn2d_dta(jbdy) == 2 ) THEN 737 if(lwp) write(numout,*) '++++++ dta%ssh/u2d/u3d allocated space' 738 IF( dta%ll_ssh ) ALLOCATE( dta%ssh(nblen(1)) ) 739 IF( dta%ll_u2d ) ALLOCATE( dta%u2d(nblen(2)) ) 740 IF( dta%ll_v2d ) ALLOCATE( dta%v2d(nblen(3)) ) 741 ENDIF 742 IF ( nn_dyn2d_dta(jbdy) == 1 .or. nn_dyn2d_dta(jbdy) == 3 ) THEN 743 IF( dta%ll_ssh ) THEN 744 if(lwp) write(numout,*) '++++++ dta%ssh pointing to fnow' 745 jfld = jfld + 1 746 dta%ssh => bf(jfld)%fnow(:,1,1) 747 ENDIF 748 IF ( dta%ll_u2d ) THEN 749 IF ( ln_full_vel_array(jbdy) ) THEN 750 if(lwp) write(numout,*) '++++++ dta%u2d allocated space' 751 ALLOCATE( dta%u2d(nblen(2)) ) 752 ELSE 753 if(lwp) write(numout,*) '++++++ dta%u2d pointing to fnow' 754 jfld = jfld + 1 755 dta%u2d => bf(jfld)%fnow(:,1,1) 756 ENDIF 757 ENDIF 758 IF ( dta%ll_v2d ) THEN 759 IF ( ln_full_vel_array(jbdy) ) THEN 760 if(lwp) write(numout,*) '++++++ dta%v2d allocated space' 761 ALLOCATE( dta%v2d(nblen(3)) ) 762 ELSE 763 if(lwp) write(numout,*) '++++++ dta%v2d pointing to fnow' 764 jfld = jfld + 1 765 dta%v2d => bf(jfld)%fnow(:,1,1) 766 ENDIF 767 ENDIF 768 ENDIF 769 770 IF ( nn_dyn3d_dta(jbdy) == 0 ) THEN 771 if(lwp) write(numout,*) '++++++ dta%u3d/v3d allocated space' 772 IF( dta%ll_u3d ) ALLOCATE( dta_bdy(jbdy)%u3d(nblen(2),jpk) ) 773 IF( dta%ll_v3d ) ALLOCATE( dta_bdy(jbdy)%v3d(nblen(3),jpk) ) 774 ENDIF 775 IF ( nn_dyn3d_dta(jbdy) == 1 .or. & 776 & ( ln_full_vel_array(jbdy) .and. ( nn_dyn2d_dta(jbdy) == 1 .or. nn_dyn2d_dta(jbdy) == 3 ) ) ) THEN 777 IF ( dta%ll_u3d .or. ( ln_full_vel_array(jbdy) .and. dta%ll_u2d ) ) THEN 778 if(lwp) write(numout,*) '++++++ dta%u3d pointing to fnow' 779 jfld = jfld + 1 780 dta_bdy(jbdy)%u3d => bf(jfld)%fnow(:,1,:) 781 ENDIF 782 IF ( dta%ll_v3d .or. ( ln_full_vel_array(jbdy) .and. dta%ll_v2d ) ) THEN 783 if(lwp) write(numout,*) '++++++ dta%v3d pointing to fnow' 784 jfld = jfld + 1 785 dta_bdy(jbdy)%v3d => bf(jfld)%fnow(:,1,:) 786 ENDIF 787 ENDIF 788 789 IF( nn_tra_dta(jbdy) == 0 ) THEN 790 if(lwp) write(numout,*) '++++++ dta%tem/sal allocated space' 791 IF( dta%ll_tem ) ALLOCATE( dta_bdy(jbdy)%tem(nblen(1),jpk) ) 792 IF( dta%ll_sal ) ALLOCATE( dta_bdy(jbdy)%sal(nblen(1),jpk) ) 793 ELSE 794 IF( dta%ll_tem ) THEN 795 if(lwp) write(numout,*) '++++++ dta%tem pointing to fnow' 796 jfld = jfld + 1 797 dta_bdy(jbdy)%tem => bf(jfld)%fnow(:,1,:) 798 ENDIF 799 IF( dta%ll_sal ) THEN 800 if(lwp) write(numout,*) '++++++ dta%sal pointing to fnow' 801 jfld = jfld + 1 802 dta_bdy(jbdy)%sal => bf(jfld)%fnow(:,1,:) 803 ENDIF 804 ENDIF 805 806 #if defined key_si3 807 IF (cn_ice(jbdy) /= 'none') THEN 808 IF( nn_ice_dta(jbdy) == 0 ) THEN 809 ALLOCATE( dta_bdy(jbdy)%a_i(nblen(1),jpl) ) 810 ALLOCATE( dta_bdy(jbdy)%h_i(nblen(1),jpl) ) 811 ALLOCATE( dta_bdy(jbdy)%h_s(nblen(1),jpl) ) 812 ELSE 813 IF ( nice_cat == jpl ) THEN ! case input cat = jpl 814 jfld = jfld + 1 815 dta_bdy(jbdy)%a_i => bf(jfld)%fnow(:,1,:) 816 jfld = jfld + 1 817 dta_bdy(jbdy)%h_i => bf(jfld)%fnow(:,1,:) 818 jfld = jfld + 1 819 dta_bdy(jbdy)%h_s => bf(jfld)%fnow(:,1,:) 820 ELSE ! case input cat = 1 OR (/=1 and /=jpl) 821 jfld_ait(jbdy) = jfld + 1 822 jfld_htit(jbdy) = jfld + 2 823 jfld_htst(jbdy) = jfld + 3 824 jfld = jfld + 3 825 ALLOCATE( dta_bdy(jbdy)%a_i(nblen(1),jpl) ) 826 ALLOCATE( dta_bdy(jbdy)%h_i(nblen(1),jpl) ) 827 ALLOCATE( dta_bdy(jbdy)%h_s(nblen(1),jpl) ) 828 dta_bdy(jbdy)%a_i(:,:) = 0._wp 829 dta_bdy(jbdy)%h_i(:,:) = 0._wp 830 dta_bdy(jbdy)%h_s(:,:) = 0._wp 831 ENDIF 832 833 ENDIF 834 ENDIF 835 #endif 517 IF( llread ) THEN ! get data from NetCDF file 518 CALL fld_fill( bf_alias, bn_alias, cn_dir, 'bdy_dta', cl3//' '//ctmp1, ctmp2 ) ! use namelist info 519 IF( bf_alias(1)%ln_tint ) ALLOCATE( bf_alias(1)%fdta( iszdim, 1, ipk, 2 ) ) 520 bf_alias(1)%imap => idx_bdy(jbdy)%nbmap(1:iszdim,igrd) ! associate the mapping used for this bdy 521 bf_alias(1)%igrd = igrd ! used only for vertical integration of 3D arrays 522 bf_alias(1)%ltotvel = ln_full_vel ! T if u3d is full velocity 523 ENDIF 524 525 ! associate the pointer and get rid of the dimensions with a size equal to 1 526 IF( jfld == jp_bdyssh ) dta_bdy(jbdy)%ssh => bf_alias(1)%fnow(:,1,1) 527 IF( jfld == jp_bdyu2d ) dta_bdy(jbdy)%u2d => bf_alias(1)%fnow(:,1,1) 528 IF( jfld == jp_bdyv2d ) dta_bdy(jbdy)%v2d => bf_alias(1)%fnow(:,1,1) 529 IF( jfld == jp_bdyu3d ) dta_bdy(jbdy)%u3d => bf_alias(1)%fnow(:,1,:) 530 IF( jfld == jp_bdyv3d ) dta_bdy(jbdy)%v3d => bf_alias(1)%fnow(:,1,:) 531 IF( jfld == jp_bdytem ) dta_bdy(jbdy)%tem => bf_alias(1)%fnow(:,1,:) 532 IF( jfld == jp_bdysal ) dta_bdy(jbdy)%sal => bf_alias(1)%fnow(:,1,:) 533 IF( jfld == jp_bdya_i ) THEN 534 IF( ipk == jpl ) THEN ; dta_bdy(jbdy)%a_i => bf_alias(1)%fnow(:,1,:) 535 ELSE ; ALLOCATE( dta_bdy(jbdy)%a_i(iszdim,jpl) ) 536 ENDIF 537 ENDIF 538 IF( jfld == jp_bdyh_i ) THEN 539 IF( ipk == jpl ) THEN ; dta_bdy(jbdy)%h_i => bf_alias(1)%fnow(:,1,:) 540 ELSE ; ALLOCATE( dta_bdy(jbdy)%h_i(iszdim,jpl) ) 541 ENDIF 542 ENDIF 543 IF( jfld == jp_bdyh_s ) THEN 544 IF( ipk == jpl ) THEN ; dta_bdy(jbdy)%h_s => bf_alias(1)%fnow(:,1,:) 545 ELSE ; ALLOCATE( dta_bdy(jbdy)%h_s(iszdim,jpl) ) 546 ENDIF 547 ENDIF 548 ENDIF 549 550 END DO ! jpbdyfld 836 551 ! 837 552 END DO ! jbdy 838 553 ! 839 554 END SUBROUTINE bdy_dta_init 840 555 841 556 !!============================================================================== 842 557 END MODULE bdydta -
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdyini.F90
r11191 r11223 37 37 38 38 INTEGER, PARAMETER :: jp_nseg = 100 ! 39 INTEGER, PARAMETER :: nrimmax = 20 ! maximum rimwidth in structured40 ! open boundary data files41 39 ! Straight open boundary segment parameters: 42 40 INTEGER :: nbdysege, nbdysegw, nbdysegn, nbdysegs … … 70 68 & cn_ice, nn_ice_dta, & 71 69 & rn_ice_tem, rn_ice_sal, rn_ice_age, & 72 & ln_vol, nn_volctl, nn_rimwidth , nb_jpk_bdy70 & ln_vol, nn_volctl, nn_rimwidth 73 71 ! 74 72 INTEGER :: ios ! Local integer output status for namelist read … … 81 79 READ ( numnam_ref, nambdy, IOSTAT = ios, ERR = 901) 82 80 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy in reference namelist', lwp ) 81 ! make sur that all elements of the namelist variables have a default definition from namelist_ref 82 ln_coords_file (2:jp_bdy) = ln_coords_file (1) 83 cn_coords_file (2:jp_bdy) = cn_coords_file (1) 84 cn_dyn2d (2:jp_bdy) = cn_dyn2d (1) 85 nn_dyn2d_dta (2:jp_bdy) = nn_dyn2d_dta (1) 86 cn_dyn3d (2:jp_bdy) = cn_dyn3d (1) 87 nn_dyn3d_dta (2:jp_bdy) = nn_dyn3d_dta (1) 88 cn_tra (2:jp_bdy) = cn_tra (1) 89 nn_tra_dta (2:jp_bdy) = nn_tra_dta (1) 90 ln_tra_dmp (2:jp_bdy) = ln_tra_dmp (1) 91 ln_dyn3d_dmp (2:jp_bdy) = ln_dyn3d_dmp (1) 92 rn_time_dmp (2:jp_bdy) = rn_time_dmp (1) 93 rn_time_dmp_out(2:jp_bdy) = rn_time_dmp_out(1) 94 cn_ice (2:jp_bdy) = cn_ice (1) 95 nn_ice_dta (2:jp_bdy) = nn_ice_dta (1) 96 rn_ice_tem (2:jp_bdy) = rn_ice_tem (1) 97 rn_ice_sal (2:jp_bdy) = rn_ice_sal (1) 98 rn_ice_age (2:jp_bdy) = rn_ice_age (1) 83 99 REWIND( numnam_cfg ) ! Namelist nambdy in configuration namelist :Unstructured open boundaries 84 100 READ ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 902 ) … … 87 103 88 104 IF( .NOT. Agrif_Root() ) ln_bdy = .FALSE. ! forced for Agrif children 105 106 IF( nb_bdy == 0 ) ln_bdy = .FALSE. 89 107 90 108 ! ----------------------------------------- … … 97 115 ! 98 116 ! Open boundaries definition (arrays and masks) 99 CALL bdy_ segs117 CALL bdy_def 100 118 ! 101 119 ! Open boundaries initialisation of external data arrays … … 115 133 116 134 117 SUBROUTINE bdy_ segs135 SUBROUTINE bdy_def 118 136 !!---------------------------------------------------------------------- 119 137 !! *** ROUTINE bdy_init *** … … 130 148 INTEGER :: ilen1, ibm1 ! - - 131 149 INTEGER :: iwe, ies, iso, ino, inum, id_dummy ! - - 132 INTEGER :: jpbdta , jpbdtau, jpbdtas! - -150 INTEGER :: jpbdta ! - - 133 151 INTEGER :: ib_bdy1, ib_bdy2, ib1, ib2 ! - - 134 152 INTEGER :: ii1, ii2, ii3, ij1, ij2, ij3 ! - - 135 153 INTEGER :: iibe, ijbe, iibi, ijbi ! - - 136 154 INTEGER :: flagu, flagv ! short cuts 155 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zz_read ! work space for 2D global boundary data 137 156 REAL(wp), POINTER, DIMENSION(:,:) :: pmask ! pointer to 2D mask fields 138 REAL(wp) :: zefl, zwfl, znfl, zsfl ! local scalars 139 INTEGER, DIMENSION (2) :: kdimsz 157 INTEGER, DIMENSION (4) :: kdimsz 140 158 INTEGER, DIMENSION(jpbgrd,jp_bdy) :: nblendta ! Length of index arrays 141 159 INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: nbidta, nbjdta ! Index arrays: i and j indices of bdy dta 142 160 INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: nbrdta ! Discrete distance from rim points 143 161 CHARACTER(LEN=1),DIMENSION(jpbgrd) :: cgrid 144 INTEGER :: jpk_max145 162 REAL(wp), DIMENSION(jpi,jpj) :: zfmask ! temporary fmask array excluding coastal boundary condition (shlat) 146 163 REAL(wp), DIMENSION(jpi,jpj) :: ztmask, zumask, zvmask ! temporary u/v mask array 147 164 !! 148 CHARACTER(LEN=1) :: ctypebdy ! - -149 165 INTEGER :: nbdyind, nbdybeg, nbdyend 150 !!151 NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend152 INTEGER :: ios ! Local integer output status for namelist read153 166 !!---------------------------------------------------------------------- 154 167 ! … … 161 174 & ' and general open boundary condition are not compatible' ) 162 175 163 IF( nb_bdy == 0 ) THEN 164 IF(lwp) WRITE(numout,*) 'nb_bdy = 0, NO OPEN BOUNDARIES APPLIED.' 165 ELSE 166 IF(lwp) WRITE(numout,*) 'Number of open boundary sets : ', nb_bdy 176 IF(lwp) WRITE(numout,*) 'Number of open boundary sets : ', nb_bdy 177 178 DO ib_bdy = 1,nb_bdy 179 180 IF(lwp) THEN 181 WRITE(numout,*) ' ' 182 WRITE(numout,*) '------ Open boundary data set ',ib_bdy,' ------' 183 IF( ln_coords_file(ib_bdy) ) THEN 184 WRITE(numout,*) 'Boundary definition read from file '//TRIM(cn_coords_file(ib_bdy)) 185 ELSE 186 WRITE(numout,*) 'Boundary defined in namelist.' 187 ENDIF 188 WRITE(numout,*) 189 ENDIF 190 191 ! barotropic bdy 192 !---------------- 193 IF(lwp) THEN 194 WRITE(numout,*) 'Boundary conditions for barotropic solution: ' 195 SELECT CASE( cn_dyn2d(ib_bdy) ) 196 CASE( 'none' ) ; WRITE(numout,*) ' no open boundary condition' 197 CASE( 'frs' ) ; WRITE(numout,*) ' Flow Relaxation Scheme' 198 CASE( 'flather' ) ; WRITE(numout,*) ' Flather radiation condition' 199 CASE( 'orlanski' ) ; WRITE(numout,*) ' Orlanski (fully oblique) radiation condition with adaptive nudging' 200 CASE( 'orlanski_npo' ) ; WRITE(numout,*) ' Orlanski (NPO) radiation condition with adaptive nudging' 201 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for cn_dyn2d' ) 202 END SELECT 203 ENDIF 204 205 dta_bdy(ib_bdy)%lneed_ssh = cn_dyn2d(ib_bdy) == 'flather' 206 dta_bdy(ib_bdy)%lneed_dyn2d = cn_dyn2d(ib_bdy) /= 'none' 207 208 IF( lwp .AND. dta_bdy(ib_bdy)%lneed_dyn2d ) THEN 209 SELECT CASE( nn_dyn2d_dta(ib_bdy) ) ! 210 CASE( 0 ) ; WRITE(numout,*) ' initial state used for bdy data' 211 CASE( 1 ) ; WRITE(numout,*) ' boundary data taken from file' 212 CASE( 2 ) ; WRITE(numout,*) ' tidal harmonic forcing taken from file' 213 CASE( 3 ) ; WRITE(numout,*) ' boundary data AND tidal harmonic forcing taken from files' 214 CASE DEFAULT ; CALL ctl_stop( 'nn_dyn2d_dta must be between 0 and 3' ) 215 END SELECT 216 ENDIF 217 IF ( dta_bdy(ib_bdy)%lneed_dyn2d .AND. nn_dyn2d_dta(ib_bdy) .GE. 2 .AND. .NOT.ln_tide ) THEN 218 CALL ctl_stop( 'You must activate with ln_tide to add tidal forcing at open boundaries' ) 219 ENDIF 220 IF(lwp) WRITE(numout,*) 221 222 ! baroclinic bdy 223 !---------------- 224 IF(lwp) THEN 225 WRITE(numout,*) 'Boundary conditions for baroclinic velocities: ' 226 SELECT CASE( cn_dyn3d(ib_bdy) ) 227 CASE('none') ; WRITE(numout,*) ' no open boundary condition' 228 CASE('frs') ; WRITE(numout,*) ' Flow Relaxation Scheme' 229 CASE('specified') ; WRITE(numout,*) ' Specified value' 230 CASE('neumann') ; WRITE(numout,*) ' Neumann conditions' 231 CASE('zerograd') ; WRITE(numout,*) ' Zero gradient for baroclinic velocities' 232 CASE('zero') ; WRITE(numout,*) ' Zero baroclinic velocities (runoff case)' 233 CASE('orlanski') ; WRITE(numout,*) ' Orlanski (fully oblique) radiation condition with adaptive nudging' 234 CASE('orlanski_npo') ; WRITE(numout,*) ' Orlanski (NPO) radiation condition with adaptive nudging' 235 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for cn_dyn3d' ) 236 END SELECT 237 ENDIF 238 239 dta_bdy(ib_bdy)%lneed_dyn3d = cn_dyn3d(ib_bdy) == 'frs' .OR. cn_dyn3d(ib_bdy) == 'specified' & 240 & .OR. cn_dyn3d(ib_bdy) == 'orlanski' .OR. cn_dyn3d(ib_bdy) == 'orlanski_npo' 241 242 IF( lwp .AND. dta_bdy(ib_bdy)%lneed_dyn3d ) THEN 243 SELECT CASE( nn_dyn3d_dta(ib_bdy) ) ! 244 CASE( 0 ) ; WRITE(numout,*) ' initial state used for bdy data' 245 CASE( 1 ) ; WRITE(numout,*) ' boundary data taken from file' 246 CASE DEFAULT ; CALL ctl_stop( 'nn_dyn3d_dta must be 0 or 1' ) 247 END SELECT 248 END IF 249 250 IF ( ln_dyn3d_dmp(ib_bdy) ) THEN 251 IF ( cn_dyn3d(ib_bdy) == 'none' ) THEN 252 IF(lwp) WRITE(numout,*) 'No open boundary condition for baroclinic velocities: ln_dyn3d_dmp is set to .false.' 253 ln_dyn3d_dmp(ib_bdy) = .false. 254 ELSEIF ( cn_dyn3d(ib_bdy) == 'frs' ) THEN 255 CALL ctl_stop( 'Use FRS OR relaxation' ) 256 ELSE 257 IF(lwp) WRITE(numout,*) ' + baroclinic velocities relaxation zone' 258 IF(lwp) WRITE(numout,*) ' Damping time scale: ',rn_time_dmp(ib_bdy),' days' 259 IF(rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' ) 260 dta_bdy(ib_bdy)%lneed_dyn3d = .TRUE. 261 ENDIF 262 ELSE 263 IF(lwp) WRITE(numout,*) ' NO relaxation on baroclinic velocities' 264 ENDIF 265 IF(lwp) WRITE(numout,*) 266 267 ! tra bdy 268 !---------------- 269 IF(lwp) THEN 270 WRITE(numout,*) 'Boundary conditions for temperature and salinity: ' 271 SELECT CASE( cn_tra(ib_bdy) ) 272 CASE('none') ; WRITE(numout,*) ' no open boundary condition' 273 CASE('frs') ; WRITE(numout,*) ' Flow Relaxation Scheme' 274 CASE('specified') ; WRITE(numout,*) ' Specified value' 275 CASE('neumann') ; WRITE(numout,*) ' Neumann conditions' 276 CASE('runoff') ; WRITE(numout,*) ' Runoff conditions : Neumann for T and specified to 0.1 for salinity' 277 CASE('orlanski') ; WRITE(numout,*) ' Orlanski (fully oblique) radiation condition with adaptive nudging' 278 CASE('orlanski_npo') ; WRITE(numout,*) ' Orlanski (NPO) radiation condition with adaptive nudging' 279 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for cn_tra' ) 280 END SELECT 281 ENDIF 282 283 dta_bdy(ib_bdy)%lneed_tra = cn_tra(ib_bdy) == 'frs' .OR. cn_tra(ib_bdy) == 'specified' & 284 & .OR. cn_tra(ib_bdy) == 'orlanski' .OR. cn_tra(ib_bdy) == 'orlanski_npo' 285 286 IF( lwp .AND. dta_bdy(ib_bdy)%lneed_tra ) THEN 287 SELECT CASE( nn_tra_dta(ib_bdy) ) ! 288 CASE( 0 ) ; WRITE(numout,*) ' initial state used for bdy data' 289 CASE( 1 ) ; WRITE(numout,*) ' boundary data taken from file' 290 CASE DEFAULT ; CALL ctl_stop( 'nn_tra_dta must be 0 or 1' ) 291 END SELECT 292 ENDIF 293 294 IF ( ln_tra_dmp(ib_bdy) ) THEN 295 IF ( cn_tra(ib_bdy) == 'none' ) THEN 296 IF(lwp) WRITE(numout,*) 'No open boundary condition for tracers: ln_tra_dmp is set to .false.' 297 ln_tra_dmp(ib_bdy) = .false. 298 ELSEIF ( cn_tra(ib_bdy) == 'frs' ) THEN 299 CALL ctl_stop( 'Use FRS OR relaxation' ) 300 ELSE 301 IF(lwp) WRITE(numout,*) ' + T/S relaxation zone' 302 IF(lwp) WRITE(numout,*) ' Damping time scale: ',rn_time_dmp(ib_bdy),' days' 303 IF(lwp) WRITE(numout,*) ' Outflow damping time scale: ',rn_time_dmp_out(ib_bdy),' days' 304 IF(lwp.AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' ) 305 dta_bdy(ib_bdy)%lneed_tra = .TRUE. 306 ENDIF 307 ELSE 308 IF(lwp) WRITE(numout,*) ' NO T/S relaxation' 309 ENDIF 310 IF(lwp) WRITE(numout,*) 311 312 #if defined key_si3 313 IF(lwp) THEN 314 WRITE(numout,*) 'Boundary conditions for sea ice: ' 315 SELECT CASE( cn_ice(ib_bdy) ) 316 CASE('none') ; WRITE(numout,*) ' no open boundary condition' 317 CASE('frs') ; WRITE(numout,*) ' Flow Relaxation Scheme' 318 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for cn_ice' ) 319 END SELECT 320 ENDIF 321 322 dta_bdy(ib_bdy)%lneed_ice = cn_ice(ib_bdy) /= 'none' 323 324 IF( lwp .AND. dta_bdy(ib_bdy)%lneed_ice ) THEN 325 SELECT CASE( nn_ice_dta(ib_bdy) ) ! 326 CASE( 0 ) ; WRITE(numout,*) ' initial state used for bdy data' 327 CASE( 1 ) ; WRITE(numout,*) ' boundary data taken from file' 328 CASE DEFAULT ; CALL ctl_stop( 'nn_ice_dta must be 0 or 1' ) 329 END SELECT 330 WRITE(numout,*) 331 WRITE(numout,*) ' tem of bdy sea-ice = ', rn_ice_tem(ib_bdy) 332 WRITE(numout,*) ' sal of bdy sea-ice = ', rn_ice_sal(ib_bdy) 333 WRITE(numout,*) ' age of bdy sea-ice = ', rn_ice_age(ib_bdy) 334 ENDIF 335 #else 336 dta_bdy(ib_bdy)%lneed_ice = .FALSE. 337 #endif 338 ! 339 IF(lwp) WRITE(numout,*) 340 IF(lwp) WRITE(numout,*) ' Width of relaxation zone = ', nn_rimwidth(ib_bdy) 341 IF(lwp) WRITE(numout,*) 342 ! 343 END DO ! nb_bdy 344 345 IF( lwp ) THEN 346 IF( ln_vol ) THEN ! check volume conservation (nn_volctl value) 347 WRITE(numout,*) 'Volume correction applied at open boundaries' 348 WRITE(numout,*) 349 SELECT CASE ( nn_volctl ) 350 CASE( 1 ) ; WRITE(numout,*) ' The total volume will be constant' 351 CASE( 0 ) ; WRITE(numout,*) ' The total volume will vary according to the surface E-P flux' 352 CASE DEFAULT ; CALL ctl_stop( 'nn_volctl must be 0 or 1' ) 353 END SELECT 354 WRITE(numout,*) 355 ! 356 ! sanity check if used with tides 357 IF( ln_tide ) THEN 358 WRITE(numout,*) ' The total volume correction is not working with tides. ' 359 WRITE(numout,*) ' Set ln_vol to .FALSE. ' 360 WRITE(numout,*) ' or ' 361 WRITE(numout,*) ' equilibriate your bdy input files ' 362 CALL ctl_stop( 'The total volume correction is not working with tides.' ) 363 END IF 364 ELSE 365 WRITE(numout,*) 'No volume correction applied at open boundaries' 366 WRITE(numout,*) 367 ENDIF 167 368 ENDIF 168 169 DO ib_bdy = 1,nb_bdy170 IF(lwp) WRITE(numout,*) ' '171 IF(lwp) WRITE(numout,*) '------ Open boundary data set ',ib_bdy,'------'172 173 IF( ln_coords_file(ib_bdy) ) THEN174 IF(lwp) WRITE(numout,*) 'Boundary definition read from file '//TRIM(cn_coords_file(ib_bdy))175 ELSE176 IF(lwp) WRITE(numout,*) 'Boundary defined in namelist.'177 ENDIF178 IF(lwp) WRITE(numout,*)179 180 IF(lwp) WRITE(numout,*) 'Boundary conditions for barotropic solution: '181 SELECT CASE( cn_dyn2d(ib_bdy) )182 CASE( 'none' )183 IF(lwp) WRITE(numout,*) ' no open boundary condition'184 dta_bdy(ib_bdy)%ll_ssh = .false.185 dta_bdy(ib_bdy)%ll_u2d = .false.186 dta_bdy(ib_bdy)%ll_v2d = .false.187 CASE( 'frs' )188 IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme'189 dta_bdy(ib_bdy)%ll_ssh = .false.190 dta_bdy(ib_bdy)%ll_u2d = .true.191 dta_bdy(ib_bdy)%ll_v2d = .true.192 CASE( 'flather' )193 IF(lwp) WRITE(numout,*) ' Flather radiation condition'194 dta_bdy(ib_bdy)%ll_ssh = .true.195 dta_bdy(ib_bdy)%ll_u2d = .true.196 dta_bdy(ib_bdy)%ll_v2d = .true.197 CASE( 'orlanski' )198 IF(lwp) WRITE(numout,*) ' Orlanski (fully oblique) radiation condition with adaptive nudging'199 dta_bdy(ib_bdy)%ll_ssh = .false.200 dta_bdy(ib_bdy)%ll_u2d = .true.201 dta_bdy(ib_bdy)%ll_v2d = .true.202 CASE( 'orlanski_npo' )203 IF(lwp) WRITE(numout,*) ' Orlanski (NPO) radiation condition with adaptive nudging'204 dta_bdy(ib_bdy)%ll_ssh = .false.205 dta_bdy(ib_bdy)%ll_u2d = .true.206 dta_bdy(ib_bdy)%ll_v2d = .true.207 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for cn_dyn2d' )208 END SELECT209 IF( cn_dyn2d(ib_bdy) /= 'none' ) THEN210 SELECT CASE( nn_dyn2d_dta(ib_bdy) ) !211 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' initial state used for bdy data'212 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' boundary data taken from file'213 CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' tidal harmonic forcing taken from file'214 CASE( 3 ) ; IF(lwp) WRITE(numout,*) ' boundary data AND tidal harmonic forcing taken from files'215 CASE DEFAULT ; CALL ctl_stop( 'nn_dyn2d_dta must be between 0 and 3' )216 END SELECT217 IF (( nn_dyn2d_dta(ib_bdy) .ge. 2 ).AND.(.NOT.ln_tide)) THEN218 CALL ctl_stop( 'You must activate with ln_tide to add tidal forcing at open boundaries' )219 ENDIF220 ENDIF221 IF(lwp) WRITE(numout,*)222 223 IF(lwp) WRITE(numout,*) 'Boundary conditions for baroclinic velocities: '224 SELECT CASE( cn_dyn3d(ib_bdy) )225 CASE('none')226 IF(lwp) WRITE(numout,*) ' no open boundary condition'227 dta_bdy(ib_bdy)%ll_u3d = .false.228 dta_bdy(ib_bdy)%ll_v3d = .false.229 CASE('frs')230 IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme'231 dta_bdy(ib_bdy)%ll_u3d = .true.232 dta_bdy(ib_bdy)%ll_v3d = .true.233 CASE('specified')234 IF(lwp) WRITE(numout,*) ' Specified value'235 dta_bdy(ib_bdy)%ll_u3d = .true.236 dta_bdy(ib_bdy)%ll_v3d = .true.237 CASE('neumann')238 IF(lwp) WRITE(numout,*) ' Neumann conditions'239 dta_bdy(ib_bdy)%ll_u3d = .false.240 dta_bdy(ib_bdy)%ll_v3d = .false.241 CASE('zerograd')242 IF(lwp) WRITE(numout,*) ' Zero gradient for baroclinic velocities'243 dta_bdy(ib_bdy)%ll_u3d = .false.244 dta_bdy(ib_bdy)%ll_v3d = .false.245 CASE('zero')246 IF(lwp) WRITE(numout,*) ' Zero baroclinic velocities (runoff case)'247 dta_bdy(ib_bdy)%ll_u3d = .false.248 dta_bdy(ib_bdy)%ll_v3d = .false.249 CASE('orlanski')250 IF(lwp) WRITE(numout,*) ' Orlanski (fully oblique) radiation condition with adaptive nudging'251 dta_bdy(ib_bdy)%ll_u3d = .true.252 dta_bdy(ib_bdy)%ll_v3d = .true.253 CASE('orlanski_npo')254 IF(lwp) WRITE(numout,*) ' Orlanski (NPO) radiation condition with adaptive nudging'255 dta_bdy(ib_bdy)%ll_u3d = .true.256 dta_bdy(ib_bdy)%ll_v3d = .true.257 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for cn_dyn3d' )258 END SELECT259 IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN260 SELECT CASE( nn_dyn3d_dta(ib_bdy) ) !261 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' initial state used for bdy data'262 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' boundary data taken from file'263 CASE DEFAULT ; CALL ctl_stop( 'nn_dyn3d_dta must be 0 or 1' )264 END SELECT265 ENDIF266 267 IF ( ln_dyn3d_dmp(ib_bdy) ) THEN268 IF ( cn_dyn3d(ib_bdy) == 'none' ) THEN269 IF(lwp) WRITE(numout,*) 'No open boundary condition for baroclinic velocities: ln_dyn3d_dmp is set to .false.'270 ln_dyn3d_dmp(ib_bdy)=.false.271 ELSEIF ( cn_dyn3d(ib_bdy) == 'frs' ) THEN272 CALL ctl_stop( 'Use FRS OR relaxation' )273 ELSE274 IF(lwp) WRITE(numout,*) ' + baroclinic velocities relaxation zone'275 IF(lwp) WRITE(numout,*) ' Damping time scale: ',rn_time_dmp(ib_bdy),' days'276 IF((lwp).AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' )277 dta_bdy(ib_bdy)%ll_u3d = .true.278 dta_bdy(ib_bdy)%ll_v3d = .true.279 ENDIF280 ELSE281 IF(lwp) WRITE(numout,*) ' NO relaxation on baroclinic velocities'282 ENDIF283 IF(lwp) WRITE(numout,*)284 285 IF(lwp) WRITE(numout,*) 'Boundary conditions for temperature and salinity: '286 SELECT CASE( cn_tra(ib_bdy) )287 CASE('none')288 IF(lwp) WRITE(numout,*) ' no open boundary condition'289 dta_bdy(ib_bdy)%ll_tem = .false.290 dta_bdy(ib_bdy)%ll_sal = .false.291 CASE('frs')292 IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme'293 dta_bdy(ib_bdy)%ll_tem = .true.294 dta_bdy(ib_bdy)%ll_sal = .true.295 CASE('specified')296 IF(lwp) WRITE(numout,*) ' Specified value'297 dta_bdy(ib_bdy)%ll_tem = .true.298 dta_bdy(ib_bdy)%ll_sal = .true.299 CASE('neumann')300 IF(lwp) WRITE(numout,*) ' Neumann conditions'301 dta_bdy(ib_bdy)%ll_tem = .false.302 dta_bdy(ib_bdy)%ll_sal = .false.303 CASE('runoff')304 IF(lwp) WRITE(numout,*) ' Runoff conditions : Neumann for T and specified to 0.1 for salinity'305 dta_bdy(ib_bdy)%ll_tem = .false.306 dta_bdy(ib_bdy)%ll_sal = .false.307 CASE('orlanski')308 IF(lwp) WRITE(numout,*) ' Orlanski (fully oblique) radiation condition with adaptive nudging'309 dta_bdy(ib_bdy)%ll_tem = .true.310 dta_bdy(ib_bdy)%ll_sal = .true.311 CASE('orlanski_npo')312 IF(lwp) WRITE(numout,*) ' Orlanski (NPO) radiation condition with adaptive nudging'313 dta_bdy(ib_bdy)%ll_tem = .true.314 dta_bdy(ib_bdy)%ll_sal = .true.315 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for cn_tra' )316 END SELECT317 IF( cn_tra(ib_bdy) /= 'none' ) THEN318 SELECT CASE( nn_tra_dta(ib_bdy) ) !319 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' initial state used for bdy data'320 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' boundary data taken from file'321 CASE DEFAULT ; CALL ctl_stop( 'nn_tra_dta must be 0 or 1' )322 END SELECT323 ENDIF324 325 IF ( ln_tra_dmp(ib_bdy) ) THEN326 IF ( cn_tra(ib_bdy) == 'none' ) THEN327 IF(lwp) WRITE(numout,*) 'No open boundary condition for tracers: ln_tra_dmp is set to .false.'328 ln_tra_dmp(ib_bdy)=.false.329 ELSEIF ( cn_tra(ib_bdy) == 'frs' ) THEN330 CALL ctl_stop( 'Use FRS OR relaxation' )331 ELSE332 IF(lwp) WRITE(numout,*) ' + T/S relaxation zone'333 IF(lwp) WRITE(numout,*) ' Damping time scale: ',rn_time_dmp(ib_bdy),' days'334 IF(lwp) WRITE(numout,*) ' Outflow damping time scale: ',rn_time_dmp_out(ib_bdy),' days'335 IF((lwp).AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' )336 dta_bdy(ib_bdy)%ll_tem = .true.337 dta_bdy(ib_bdy)%ll_sal = .true.338 ENDIF339 ELSE340 IF(lwp) WRITE(numout,*) ' NO T/S relaxation'341 ENDIF342 IF(lwp) WRITE(numout,*)343 344 #if defined key_si3345 IF(lwp) WRITE(numout,*) 'Boundary conditions for sea ice: '346 SELECT CASE( cn_ice(ib_bdy) )347 CASE('none')348 IF(lwp) WRITE(numout,*) ' no open boundary condition'349 dta_bdy(ib_bdy)%ll_a_i = .false.350 dta_bdy(ib_bdy)%ll_h_i = .false.351 dta_bdy(ib_bdy)%ll_h_s = .false.352 CASE('frs')353 IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme'354 dta_bdy(ib_bdy)%ll_a_i = .true.355 dta_bdy(ib_bdy)%ll_h_i = .true.356 dta_bdy(ib_bdy)%ll_h_s = .true.357 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for cn_ice' )358 END SELECT359 IF( cn_ice(ib_bdy) /= 'none' ) THEN360 SELECT CASE( nn_ice_dta(ib_bdy) ) !361 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' initial state used for bdy data'362 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' boundary data taken from file'363 CASE DEFAULT ; CALL ctl_stop( 'nn_ice_dta must be 0 or 1' )364 END SELECT365 ENDIF366 IF(lwp) WRITE(numout,*)367 IF(lwp) WRITE(numout,*) ' tem of bdy sea-ice = ', rn_ice_tem(ib_bdy)368 IF(lwp) WRITE(numout,*) ' sal of bdy sea-ice = ', rn_ice_sal(ib_bdy)369 IF(lwp) WRITE(numout,*) ' age of bdy sea-ice = ', rn_ice_age(ib_bdy)370 #endif371 372 IF(lwp) WRITE(numout,*) ' Width of relaxation zone = ', nn_rimwidth(ib_bdy)373 IF(lwp) WRITE(numout,*)374 !375 END DO376 377 IF( nb_bdy > 0 ) THEN378 IF( ln_vol ) THEN ! check volume conservation (nn_volctl value)379 IF(lwp) WRITE(numout,*) 'Volume correction applied at open boundaries'380 IF(lwp) WRITE(numout,*)381 SELECT CASE ( nn_volctl )382 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' The total volume will be constant'383 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' The total volume will vary according to the surface E-P flux'384 CASE DEFAULT ; CALL ctl_stop( 'nn_volctl must be 0 or 1' )385 END SELECT386 IF(lwp) WRITE(numout,*)387 !388 ! sanity check if used with tides389 IF( ln_tide ) THEN390 IF(lwp) WRITE(numout,*) ' The total volume correction is not working with tides. '391 IF(lwp) WRITE(numout,*) ' Set ln_vol to .FALSE. '392 IF(lwp) WRITE(numout,*) ' or '393 IF(lwp) WRITE(numout,*) ' equilibriate your bdy input files '394 CALL ctl_stop( 'The total volume correction is not working with tides.' )395 END IF396 ELSE397 IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries'398 IF(lwp) WRITE(numout,*)399 ENDIF400 IF( nb_jpk_bdy(ib_bdy) > 0 ) THEN401 IF(lwp) WRITE(numout,*) '*** open boundary will be interpolate in the vertical onto the native grid ***'402 ELSE403 IF(lwp) WRITE(numout,*) '*** open boundary will be read straight onto the native grid without vertical interpolation ***'404 ENDIF405 ENDIF406 369 407 370 ! ------------------------------------------------- … … 409 372 ! ------------------------------------------------- 410 373 411 ! Work out global dimensions of boundary data412 ! ---------------------------------------------413 374 REWIND( numnam_cfg ) 414 415 375 nblendta(:,:) = 0 416 376 nbdysege = 0 … … 418 378 nbdysegn = 0 419 379 nbdysegs = 0 420 icount = 0 ! count user defined segments 421 ! Dimensions below are used to allocate arrays to read external data 422 jpbdtas = 1 ! Maximum size of boundary data (structured case) 423 jpbdtau = 1 ! Maximum size of boundary data (unstructured case) 424 380 381 ! Define all boundaries 382 ! --------------------- 425 383 DO ib_bdy = 1, nb_bdy 426 427 IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Work out size of global arrays from namelist parameters 428 429 icount = icount + 1 430 ! No REWIND here because may need to read more than one nambdy_index namelist. 431 ! Read only namelist_cfg to avoid unseccessfull overwrite 432 ! keep full control of the configuration namelist 433 READ ( numnam_cfg, nambdy_index, IOSTAT = ios, ERR = 904 ) 434 904 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy_index in configuration namelist', lwp ) 435 IF(lwm) WRITE ( numond, nambdy_index ) 436 437 SELECT CASE ( TRIM(ctypebdy) ) 438 CASE( 'N' ) 439 IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1 440 nbdyind = jpjglo - 2 ! set boundary to whole side of model domain. 441 nbdybeg = 2 442 nbdyend = jpiglo - 1 443 ENDIF 444 nbdysegn = nbdysegn + 1 445 npckgn(nbdysegn) = ib_bdy ! Save bdy package number 446 jpjnob(nbdysegn) = nbdyind 447 jpindt(nbdysegn) = nbdybeg 448 jpinft(nbdysegn) = nbdyend 449 ! 450 CASE( 'S' ) 451 IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1 452 nbdyind = 2 ! set boundary to whole side of model domain. 453 nbdybeg = 2 454 nbdyend = jpiglo - 1 455 ENDIF 456 nbdysegs = nbdysegs + 1 457 npckgs(nbdysegs) = ib_bdy ! Save bdy package number 458 jpjsob(nbdysegs) = nbdyind 459 jpisdt(nbdysegs) = nbdybeg 460 jpisft(nbdysegs) = nbdyend 461 ! 462 CASE( 'E' ) 463 IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1 464 nbdyind = jpiglo - 2 ! set boundary to whole side of model domain. 465 nbdybeg = 2 466 nbdyend = jpjglo - 1 467 ENDIF 468 nbdysege = nbdysege + 1 469 npckge(nbdysege) = ib_bdy ! Save bdy package number 470 jpieob(nbdysege) = nbdyind 471 jpjedt(nbdysege) = nbdybeg 472 jpjeft(nbdysege) = nbdyend 473 ! 474 CASE( 'W' ) 475 IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1 476 nbdyind = 2 ! set boundary to whole side of model domain. 477 nbdybeg = 2 478 nbdyend = jpjglo - 1 479 ENDIF 480 nbdysegw = nbdysegw + 1 481 npckgw(nbdysegw) = ib_bdy ! Save bdy package number 482 jpiwob(nbdysegw) = nbdyind 483 jpjwdt(nbdysegw) = nbdybeg 484 jpjwft(nbdysegw) = nbdyend 485 ! 486 CASE DEFAULT ; CALL ctl_stop( 'ctypebdy must be N, S, E or W' ) 487 END SELECT 488 489 ! For simplicity we assume that in case of straight bdy, arrays have the same length 490 ! (even if it is true that last tangential velocity points 491 ! are useless). This simplifies a little bit boundary data format (and agrees with format 492 ! used so far in obc package) 493 494 nblendta(1:jpbgrd,ib_bdy) = (nbdyend - nbdybeg + 1) * nn_rimwidth(ib_bdy) 495 jpbdtas = MAX(jpbdtas, (nbdyend - nbdybeg + 1)) 496 IF (lwp.and.(nn_rimwidth(ib_bdy)>nrimmax)) & 497 & CALL ctl_stop( 'rimwidth must be lower than nrimmax' ) 498 499 ELSE ! Read size of arrays in boundary coordinates file. 384 ! 385 IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! build bdy coordinates with segments defined in namelist 386 387 CALL bdy_read_seg( ib_bdy, nblendta(:,ib_bdy) ) 388 ! Now look for crossings in user (namelist) defined open boundary segments: 389 IF( nbdysege > 0 .OR. nbdysegw > 0 .OR. nbdysegn > 0 .OR. nbdysegs > 0) CALL bdy_ctl_seg 390 391 ELSE ! Read size of arrays in boundary coordinates file. 392 500 393 CALL iom_open( cn_coords_file(ib_bdy), inum ) 501 394 DO igrd = 1, jpbgrd 502 395 id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz ) 503 396 nblendta(igrd,ib_bdy) = MAXVAL(kdimsz) 504 jpbdtau = MAX(jpbdtau, MAXVAL(kdimsz))505 397 END DO 506 398 CALL iom_close( inum ) 507 ! 508 ENDIF 399 ENDIF 509 400 ! 510 401 END DO ! ib_bdy 511 512 IF (nb_bdy>0) THEN 513 jpbdta = MAXVAL(nblendta(1:jpbgrd,1:nb_bdy)) 514 515 ! Allocate arrays 516 !--------------- 517 ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy), & 518 & nbrdta(jpbdta, jpbgrd, nb_bdy) ) 519 520 jpk_max = MAXVAL(nb_jpk_bdy) 521 jpk_max = MAX(jpk_max, jpk) 522 523 ALLOCATE( dta_global(jpbdtau, 1, jpk_max) ) 524 ALLOCATE( dta_global_z(jpbdtau, 1, jpk_max) ) ! needed ?? TODO 525 ALLOCATE( dta_global_dz(jpbdtau, 1, jpk_max) )! needed ?? TODO 526 527 IF ( icount>0 ) THEN 528 ALLOCATE( dta_global2(jpbdtas, nrimmax, jpk_max) ) 529 ALLOCATE( dta_global2_z(jpbdtas, nrimmax, jpk_max) ) ! needed ?? TODO 530 ALLOCATE( dta_global2_dz(jpbdtas, nrimmax, jpk_max) )! needed ?? TODO 531 ENDIF 532 ! 533 ENDIF 534 535 ! Now look for crossings in user (namelist) defined open boundary segments: 536 !-------------------------------------------------------------------------- 537 IF( icount>0 ) CALL bdy_ctl_seg 538 402 403 ! Allocate arrays 404 !--------------- 405 jpbdta = MAXVAL(nblendta(1:jpbgrd,1:nb_bdy)) 406 ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy), nbrdta(jpbdta, jpbgrd, nb_bdy) ) 407 539 408 ! Calculate global boundary index arrays or read in from file 540 409 !------------------------------------------------------------ … … 544 413 IF( ln_coords_file(ib_bdy) ) THEN 545 414 ! 415 ALLOCATE( zz_read( MAXVAL(nblendta), 1 ) ) 546 416 CALL iom_open( cn_coords_file(ib_bdy), inum ) 417 ! 547 418 DO igrd = 1, jpbgrd 548 CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )419 CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), zz_read(1:nblendta(igrd,ib_bdy),:) ) 549 420 DO ii = 1,nblendta(igrd,ib_bdy) 550 nbidta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )421 nbidta(ii,igrd,ib_bdy) = NINT( zz_read(ii,1) ) 551 422 END DO 552 CALL iom_get( inum, jpdom_unknown, 'nbj'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )423 CALL iom_get( inum, jpdom_unknown, 'nbj'//cgrid(igrd), zz_read(1:nblendta(igrd,ib_bdy),:) ) 553 424 DO ii = 1,nblendta(igrd,ib_bdy) 554 nbjdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )425 nbjdta(ii,igrd,ib_bdy) = NINT( zz_read(ii,1) ) 555 426 END DO 556 CALL iom_get( inum, jpdom_unknown, 'nbr'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )427 CALL iom_get( inum, jpdom_unknown, 'nbr'//cgrid(igrd), zz_read(1:nblendta(igrd,ib_bdy),:) ) 557 428 DO ii = 1,nblendta(igrd,ib_bdy) 558 nbrdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )429 nbrdta(ii,igrd,ib_bdy) = NINT( zz_read(ii,1) ) 559 430 END DO 560 431 ! … … 564 435 IF(lwp) WRITE(numout,*) ' nn_rimwidth from namelist is ', nn_rimwidth(ib_bdy) 565 436 IF (ibr_max < nn_rimwidth(ib_bdy)) & 566 CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_bdy) ) 567 END DO 437 CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_bdy) ) 438 END DO 439 ! 568 440 CALL iom_close( inum ) 441 DEALLOCATE( zz_read ) 569 442 ! 570 ENDIF 571 ! 572 END DO 573 443 ENDIF 444 ! 445 END DO 446 574 447 ! 2. Now fill indices corresponding to straight open boundary arrays: 575 ! East 576 !----- 577 DO iseg = 1, nbdysege 578 ib_bdy = npckge(iseg) 579 ! 580 ! ------------ T points ------------- 581 igrd=1 582 icount=0 583 DO ir = 1, nn_rimwidth(ib_bdy) 584 DO ij = jpjedt(iseg), jpjeft(iseg) 585 icount = icount + 1 586 nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir 587 nbjdta(icount, igrd, ib_bdy) = ij 588 nbrdta(icount, igrd, ib_bdy) = ir 589 ENDDO 590 ENDDO 591 ! 592 ! ------------ U points ------------- 593 igrd=2 594 icount=0 595 DO ir = 1, nn_rimwidth(ib_bdy) 596 DO ij = jpjedt(iseg), jpjeft(iseg) 597 icount = icount + 1 598 nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 1 - ir 599 nbjdta(icount, igrd, ib_bdy) = ij 600 nbrdta(icount, igrd, ib_bdy) = ir 601 ENDDO 602 ENDDO 603 ! 604 ! ------------ V points ------------- 605 igrd=3 606 icount=0 607 DO ir = 1, nn_rimwidth(ib_bdy) 608 ! DO ij = jpjedt(iseg), jpjeft(iseg) - 1 609 DO ij = jpjedt(iseg), jpjeft(iseg) 610 icount = icount + 1 611 nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir 612 nbjdta(icount, igrd, ib_bdy) = ij 613 nbrdta(icount, igrd, ib_bdy) = ir 614 ENDDO 615 nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 616 nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 617 ENDDO 618 ENDDO 619 ! 620 ! West 621 !----- 622 DO iseg = 1, nbdysegw 623 ib_bdy = npckgw(iseg) 624 ! 625 ! ------------ T points ------------- 626 igrd=1 627 icount=0 628 DO ir = 1, nn_rimwidth(ib_bdy) 629 DO ij = jpjwdt(iseg), jpjwft(iseg) 630 icount = icount + 1 631 nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 632 nbjdta(icount, igrd, ib_bdy) = ij 633 nbrdta(icount, igrd, ib_bdy) = ir 634 ENDDO 635 ENDDO 636 ! 637 ! ------------ U points ------------- 638 igrd=2 639 icount=0 640 DO ir = 1, nn_rimwidth(ib_bdy) 641 DO ij = jpjwdt(iseg), jpjwft(iseg) 642 icount = icount + 1 643 nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 644 nbjdta(icount, igrd, ib_bdy) = ij 645 nbrdta(icount, igrd, ib_bdy) = ir 646 ENDDO 647 ENDDO 648 ! 649 ! ------------ V points ------------- 650 igrd=3 651 icount=0 652 DO ir = 1, nn_rimwidth(ib_bdy) 653 ! DO ij = jpjwdt(iseg), jpjwft(iseg) - 1 654 DO ij = jpjwdt(iseg), jpjwft(iseg) 655 icount = icount + 1 656 nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 657 nbjdta(icount, igrd, ib_bdy) = ij 658 nbrdta(icount, igrd, ib_bdy) = ir 659 ENDDO 660 nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 661 nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 662 ENDDO 663 ENDDO 664 ! 665 ! North 666 !----- 667 DO iseg = 1, nbdysegn 668 ib_bdy = npckgn(iseg) 669 ! 670 ! ------------ T points ------------- 671 igrd=1 672 icount=0 673 DO ir = 1, nn_rimwidth(ib_bdy) 674 DO ii = jpindt(iseg), jpinft(iseg) 675 icount = icount + 1 676 nbidta(icount, igrd, ib_bdy) = ii 677 nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir 678 nbrdta(icount, igrd, ib_bdy) = ir 679 ENDDO 680 ENDDO 681 ! 682 ! ------------ U points ------------- 683 igrd=2 684 icount=0 685 DO ir = 1, nn_rimwidth(ib_bdy) 686 ! DO ii = jpindt(iseg), jpinft(iseg) - 1 687 DO ii = jpindt(iseg), jpinft(iseg) 688 icount = icount + 1 689 nbidta(icount, igrd, ib_bdy) = ii 690 nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir 691 nbrdta(icount, igrd, ib_bdy) = ir 692 ENDDO 693 nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 694 nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 695 ENDDO 696 ! 697 ! ------------ V points ------------- 698 igrd=3 699 icount=0 700 DO ir = 1, nn_rimwidth(ib_bdy) 701 DO ii = jpindt(iseg), jpinft(iseg) 702 icount = icount + 1 703 nbidta(icount, igrd, ib_bdy) = ii 704 nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 1 - ir 705 nbrdta(icount, igrd, ib_bdy) = ir 706 ENDDO 707 ENDDO 708 ENDDO 709 ! 710 ! South 711 !----- 712 DO iseg = 1, nbdysegs 713 ib_bdy = npckgs(iseg) 714 ! 715 ! ------------ T points ------------- 716 igrd=1 717 icount=0 718 DO ir = 1, nn_rimwidth(ib_bdy) 719 DO ii = jpisdt(iseg), jpisft(iseg) 720 icount = icount + 1 721 nbidta(icount, igrd, ib_bdy) = ii 722 nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 723 nbrdta(icount, igrd, ib_bdy) = ir 724 ENDDO 725 ENDDO 726 ! 727 ! ------------ U points ------------- 728 igrd=2 729 icount=0 730 DO ir = 1, nn_rimwidth(ib_bdy) 731 ! DO ii = jpisdt(iseg), jpisft(iseg) - 1 732 DO ii = jpisdt(iseg), jpisft(iseg) 733 icount = icount + 1 734 nbidta(icount, igrd, ib_bdy) = ii 735 nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 736 nbrdta(icount, igrd, ib_bdy) = ir 737 ENDDO 738 nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 739 nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 740 ENDDO 741 ! 742 ! ------------ V points ------------- 743 igrd=3 744 icount=0 745 DO ir = 1, nn_rimwidth(ib_bdy) 746 DO ii = jpisdt(iseg), jpisft(iseg) 747 icount = icount + 1 748 nbidta(icount, igrd, ib_bdy) = ii 749 nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 750 nbrdta(icount, igrd, ib_bdy) = ir 751 ENDDO 752 ENDDO 753 ENDDO 448 CALL bdy_coords_seg( nbidta, nbjdta, nbrdta ) 754 449 755 450 ! Deal with duplicated points … … 765 460 DO ib2 = 1, nblendta(igrd,ib_bdy2) 766 461 IF ((nbidta(ib1, igrd, ib_bdy1)==nbidta(ib2, igrd, ib_bdy2)).AND. & 767 & (nbjdta(ib1, igrd, ib_bdy1)==nbjdta(ib2, igrd, ib_bdy2))) THEN768 ! IF ((lwp).AND.(igrd==1)) WRITE(numout,*) ' found coincident point ji, jj:', &769 ! & nbidta(ib1, igrd, ib_bdy1), &770 ! & nbjdta(ib2, igrd, ib_bdy2)462 & (nbjdta(ib1, igrd, ib_bdy1)==nbjdta(ib2, igrd, ib_bdy2))) THEN 463 ! IF ((lwp).AND.(igrd==1)) WRITE(numout,*) ' found coincident point ji, jj:', & 464 ! & nbidta(ib1, igrd, ib_bdy1), & 465 ! & nbjdta(ib2, igrd, ib_bdy2) 771 466 ! keep only points with the lowest distance to boundary: 772 467 IF (nbrdta(ib1, igrd, ib_bdy1)<nbrdta(ib2, igrd, ib_bdy2)) THEN 773 nbidta(ib2, igrd, ib_bdy2) =-ib_bdy2774 nbjdta(ib2, igrd, ib_bdy2) =-ib_bdy2468 nbidta(ib2, igrd, ib_bdy2) =-ib_bdy2 469 nbjdta(ib2, igrd, ib_bdy2) =-ib_bdy2 775 470 ELSEIF (nbrdta(ib1, igrd, ib_bdy1)>nbrdta(ib2, igrd, ib_bdy2)) THEN 776 nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1777 nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1778 ! Arbitrary choice if distances are the same:471 nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1 472 nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1 473 ! Arbitrary choice if distances are the same: 779 474 ELSE 780 nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1781 nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1475 nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1 476 nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1 782 477 ENDIF 783 478 END IF … … 811 506 IF( nbrdta(ib,igrd,ib_bdy) < nbrdta(ibm1,igrd,ib_bdy) ) THEN 812 507 CALL ctl_stop('bdy_segs : ERROR : boundary data in file must be defined ', & 813 814 815 ENDIF 508 & ' in order of distance from edge nbr A utility for re-ordering ', & 509 & ' boundary coordinates and data files exists in the TOOLS/OBC directory') 510 ENDIF 816 511 ENDIF 817 512 ! check if point is in local domain … … 911 606 ir = MAX( 1, idx_bdy(ib_bdy)%nbr(ib,igrd) ) ! both rim 0 and rim 1 have the same weights 912 607 idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( REAL( ir - 1 ) *0.5 ) ! tanh formulation 913 ! idx_bdy(ib_bdy)%nbw(ib,igrd) = (REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2. ! quadratic914 ! idx_bdy(ib_bdy)%nbw(ib,igrd) = REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)) ! linear915 END DO 916 END DO 608 ! idx_bdy(ib_bdy)%nbw(ib,igrd) = (REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2. ! quadratic 609 ! idx_bdy(ib_bdy)%nbw(ib,igrd) = REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)) ! linear 610 END DO 611 END DO 917 612 918 613 ! Compute damping coefficients … … 922 617 ir = MAX( 1, idx_bdy(ib_bdy)%nbr(ib,igrd) ) ! both rim 0 and rim 1 have the same damping coefficients 923 618 idx_bdy(ib_bdy)%nbd(ib,igrd) = 1. / ( rn_time_dmp(ib_bdy) * rday ) & 924 & *(REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2. ! quadratic619 & *(REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2. ! quadratic 925 620 idx_bdy(ib_bdy)%nbdout(ib,igrd) = 1. / ( rn_time_dmp_out(ib_bdy) * rday ) & 926 & *(REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2. ! quadratic927 END DO 928 END DO 621 & *(REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2. ! quadratic 622 END DO 623 END DO 929 624 930 625 END DO ! ib_bdy … … 935 630 936 631 ! ------------------------------------------ 937 ! hand elrim0, do as if rim 1 was free ocean632 ! handle rim0, do as if rim 1 was free ocean 938 633 ! ------------------------------------------ 939 634 … … 944 639 DO ii = 2, jpim1 945 640 zfmask(ii,ij) = ztmask(ii,ij ) * ztmask(ii+1,ij ) & 946 & * ztmask(ii,ij+1) * ztmask(ii+1,ij+1)947 END DO 641 & * ztmask(ii,ij+1) * ztmask(ii+1,ij+1) 642 END DO 948 643 END DO 949 644 CALL lbc_lnk( 'bdyini', zfmask, 'F', 1. ) … … 953 648 ! bdytmask = 1 on the computational domain AND on open boundaries 954 649 ! = 0 elsewhere 955 650 956 651 bdytmask(:,:) = ssmask(:,:) 957 652 … … 981 676 982 677 ! ------------------------------------ 983 ! hand elrim1, do as if rim 0 was land678 ! handle rim1, do as if rim 0 was land 984 679 ! ------------------------------------ 985 680 … … 1000 695 DO ii = 2, jpim1 1001 696 zfmask(ii,ij) = ztmask(ii,ij ) * ztmask(ii+1,ij ) & 1002 & * ztmask(ii,ij+1) * ztmask(ii+1,ij+1)1003 END DO 697 & * ztmask(ii,ij+1) * ztmask(ii+1,ij+1) 698 END DO 1004 699 END DO 1005 700 CALL lbc_lnk( 'bdyini', zfmask, 'F', 1. ) … … 1019 714 1020 715 CALL bdy_rim_treat( zumask, zvmask, zfmask, .false. ) ! compute flagu, flagv, ntreat on rim 1 1021 1022 716 1023 717 ! … … 1053 747 IF( iibi == 0 .OR. ii1 == 0 .OR. ii2 == 0 .OR. ii3 == 0 ) lrecv_bdyint(ib_bdy,igrd,1,ir) = .true. 1054 748 IF( iibi == jpi+1 .OR. ii1 == jpi+1 .OR. ii2 == jpi+1 .OR. ii3 == jpi+1 ) lrecv_bdyint(ib_bdy,igrd,2,ir) = .true. 1055 IF( iibe == 0 ) lrecv_bdyext(ib_bdy,igrd,1,ir) = .true.1056 IF( iibe == jpi+1 ) lrecv_bdyext(ib_bdy,igrd,2,ir) = .true.749 IF( iibe == 0 ) lrecv_bdyext(ib_bdy,igrd,1,ir) = .true. 750 IF( iibe == jpi+1 ) lrecv_bdyext(ib_bdy,igrd,2,ir) = .true. 1057 751 ! Check if neighbour has its rim parallel to its mpi subdomain border and located next to its halo 1058 752 ! :¨¨¨¨¨|¨¨--> | | <--¨¨|¨¨¨¨¨: … … 1060 754 ! :.....|_._:_____| (1) W neighbour E neighbour (2) |_____:_._|.....: 1061 755 IF( ii == 2 .AND. ( nbondi == 1 .OR. nbondi == 0 ) .AND. & 1062 & ( iibi == 3 .OR. ii1 == 3 .OR. ii2 == 3 .OR. ii3 == 3 ) ) lsend_bdyint(ib_bdy,igrd,1,ir) =.true.756 & ( iibi == 3 .OR. ii1 == 3 .OR. ii2 == 3 .OR. ii3 == 3 ) ) lsend_bdyint(ib_bdy,igrd,1,ir)=.true. 1063 757 IF( ii == jpi-1 .AND. ( nbondi == -1 .OR. nbondi == 0 ) .AND. & 1064 & ( iibi == jpi-2 .OR. ii1 == jpi-2 .OR. ii2 == jpi-2 .OR. ii3 == jpi-2) ) lsend_bdyint(ib_bdy,igrd,2,ir) =.true.1065 IF( ii == 2 .AND. ( nbondi == 1 .OR. nbondi == 0 ) .AND. iibe == 3 ) lsend_bdyext(ib_bdy,igrd,1,ir) =.true.1066 IF( ii == jpi-1 .AND. ( nbondi == -1 .OR. nbondi == 0 ) .AND. iibe == jpi-2 ) lsend_bdyext(ib_bdy,igrd,2,ir) =.true.758 & ( iibi == jpi-2 .OR. ii1 == jpi-2 .OR. ii2 == jpi-2 .OR. ii3 == jpi-2) ) lsend_bdyint(ib_bdy,igrd,2,ir)=.true. 759 IF( ii == 2 .AND. ( nbondi == 1 .OR. nbondi == 0 ) .AND. iibe == 3 ) lsend_bdyext(ib_bdy,igrd,1,ir)=.true. 760 IF( ii == jpi-1 .AND. ( nbondi == -1 .OR. nbondi == 0 ) .AND. iibe == jpi-2 ) lsend_bdyext(ib_bdy,igrd,2,ir)=.true. 1067 761 ! 1068 762 ! search neighbour in the north/south direction … … 1073 767 IF( ijbi == 0 .OR. ij1 == 0 .OR. ij2 == 0 .OR. ij3 == 0 ) lrecv_bdyint(ib_bdy,igrd,3,ir) = .true. 1074 768 IF( ijbi == jpj+1 .OR. ij1 == jpj+1 .OR. ij2 == jpj+1 .OR. ij3 == jpj+1 ) lrecv_bdyint(ib_bdy,igrd,4,ir) = .true. 1075 IF( ijbe == 0 ) lrecv_bdyext(ib_bdy,igrd,3,ir) = .true.1076 IF( ijbe == jpj+1 ) lrecv_bdyext(ib_bdy,igrd,4,ir) = .true.769 IF( ijbe == 0 ) lrecv_bdyext(ib_bdy,igrd,3,ir) = .true. 770 IF( ijbe == jpj+1 ) lrecv_bdyext(ib_bdy,igrd,4,ir) = .true. 1077 771 ! Check if neighbour has its rim parallel to its mpi subdomain _________ border and next to its halo 1078 772 ! ^ | o | : : … … 1080 774 ! :_________: (3) S neighbour N neighbour (4) v | o | 1081 775 IF( ij == 2 .AND. ( nbondj == 1 .OR. nbondj == 0 ) .AND. & 1082 & ( ijbi == 3 .OR. ij1 == 3 .OR. ij2 == 3 .OR. ij3 == 3 ) ) lsend_bdyint(ib_bdy,igrd,3,ir) =.true.776 & ( ijbi == 3 .OR. ij1 == 3 .OR. ij2 == 3 .OR. ij3 == 3 ) ) lsend_bdyint(ib_bdy,igrd,3,ir)=.true. 1083 777 IF( ij == jpj-1 .AND. ( nbondj == -1 .OR. nbondj == 0 ) .AND. & 1084 & ( ijbi == jpj-2 .OR. ij1 == jpj-2 .OR. ij2 == jpj-2 .OR. ij3 == jpj-2) ) lsend_bdyint(ib_bdy,igrd,4,ir) =.true.1085 IF( ij == 2 .AND. ( nbondj == 1 .OR. nbondj == 0 ) .AND. ijbe == 3 ) lsend_bdyext(ib_bdy,igrd,3,ir) =.true.1086 IF( ij == jpj-1 .AND. ( nbondj == -1 .OR. nbondj == 0 ) .AND. ijbe == jpj-2 ) lsend_bdyext(ib_bdy,igrd,4,ir) =.true.778 & ( ijbi == jpj-2 .OR. ij1 == jpj-2 .OR. ij2 == jpj-2 .OR. ij3 == jpj-2) ) lsend_bdyint(ib_bdy,igrd,4,ir)=.true. 779 IF( ij == 2 .AND. ( nbondj == 1 .OR. nbondj == 0 ) .AND. ijbe == 3 ) lsend_bdyext(ib_bdy,igrd,3,ir)=.true. 780 IF( ij == jpj-1 .AND. ( nbondj == -1 .OR. nbondj == 0 ) .AND. ijbe == jpj-2 ) lsend_bdyext(ib_bdy,igrd,4,ir)=.true. 1087 781 END DO 1088 782 END DO … … 1091 785 DO ib_bdy = 1,nb_bdy 1092 786 IF( cn_dyn2d(ib_bdy) == 'orlanski' .OR. cn_dyn2d(ib_bdy) == 'orlanski_npo' .OR. & 1093 & cn_dyn3d(ib_bdy) == 'orlanski' .OR. cn_dyn3d(ib_bdy) == 'orlanski_npo' .OR. &1094 & cn_tra(ib_bdy) == 'orlanski' .OR. cn_tra(ib_bdy) == 'orlanski_npo' ) THEN787 & cn_dyn3d(ib_bdy) == 'orlanski' .OR. cn_dyn3d(ib_bdy) == 'orlanski_npo' .OR. & 788 & cn_tra(ib_bdy) == 'orlanski' .OR. cn_tra(ib_bdy) == 'orlanski_npo' ) THEN 1095 789 DO igrd = 1, jpbgrd 1096 790 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) … … 1110 804 ! Tidy up 1111 805 !-------- 1112 IF( nb_bdy>0 )DEALLOCATE( nbidta, nbjdta, nbrdta )1113 ! 1114 END SUBROUTINE bdy_ segs806 DEALLOCATE( nbidta, nbjdta, nbrdta ) 807 ! 808 END SUBROUTINE bdy_def 1115 809 1116 810 … … 1120 814 !! 1121 815 !! ** Purpose : Initialize structures ( flagu, flagv, ntreat ) indicating how rim points 1122 !! are to be hand eled in the boundary condition treatment1123 !! 1124 !! ** Method : - to hand elrim 0 zmasks must indicate ocean points (set at one on rim 0 and rim 1 and interior)816 !! are to be handled in the boundary condition treatment 817 !! 818 !! ** Method : - to handle rim 0 zmasks must indicate ocean points (set at one on rim 0 and rim 1 and interior) 1125 819 !! and bdymasks must be set at 0 on rim 0 (set at one on rim 1 and interior) 1126 820 !! (as if rim 1 was free ocean) 1127 !! - to hand elrim 1 zmasks must be set at 0 on rim 0 (set at one on rim 1 and interior)821 !! - to handle rim 1 zmasks must be set at 0 on rim 0 (set at one on rim 1 and interior) 1128 822 !! and bdymasks must indicate free ocean points (set at one on interior) 1129 823 !! (as if rim 0 was land) … … 1373 1067 CASE( 16 ) ; ii1 = ii-1 ; ij1 = ij ; ii2 = ii ; ij2 = ij+1 ; ii3 = ii+1 ; ij3 = ij 1374 1068 END SELECT 1375 1069 END SUBROUTINE find_neib 1376 1070 1377 1071 1072 SUBROUTINE bdy_read_seg( kb_bdy, knblendta ) 1073 !!---------------------------------------------------------------------- 1074 !! *** ROUTINE bdy_coords_seg *** 1075 !! 1076 !! ** Purpose : build bdy coordinates with segments defined in namelist 1077 !! 1078 !! ** Method : read namelist nambdy_index blocks 1079 !! 1080 !!---------------------------------------------------------------------- 1081 INTEGER , INTENT (in ) :: kb_bdy ! bdy number 1082 INTEGER, DIMENSION(jpbgrd), INTENT ( out) :: knblendta ! length of index arrays 1083 !! 1084 INTEGER :: ios ! Local integer output status for namelist read 1085 INTEGER :: nbdyind, nbdybeg, nbdyend 1086 CHARACTER(LEN=1) :: ctypebdy ! - - 1087 NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend 1088 !!---------------------------------------------------------------------- 1089 1090 ! No REWIND here because may need to read more than one nambdy_index namelist. 1091 ! Read only namelist_cfg to avoid unseccessfull overwrite 1092 ! keep full control of the configuration namelist 1093 READ ( numnam_cfg, nambdy_index, IOSTAT = ios, ERR = 904 ) 1094 904 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy_index in configuration namelist', lwp ) 1095 IF(lwm) WRITE ( numond, nambdy_index ) 1096 1097 SELECT CASE ( TRIM(ctypebdy) ) 1098 CASE( 'N' ) 1099 IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1 1100 nbdyind = jpjglo - 2 ! set boundary to whole side of model domain. 1101 nbdybeg = 2 1102 nbdyend = jpiglo - 1 1103 ENDIF 1104 nbdysegn = nbdysegn + 1 1105 npckgn(nbdysegn) = kb_bdy ! Save bdy package number 1106 jpjnob(nbdysegn) = nbdyind 1107 jpindt(nbdysegn) = nbdybeg 1108 jpinft(nbdysegn) = nbdyend 1109 ! 1110 CASE( 'S' ) 1111 IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1 1112 nbdyind = 2 ! set boundary to whole side of model domain. 1113 nbdybeg = 2 1114 nbdyend = jpiglo - 1 1115 ENDIF 1116 nbdysegs = nbdysegs + 1 1117 npckgs(nbdysegs) = kb_bdy ! Save bdy package number 1118 jpjsob(nbdysegs) = nbdyind 1119 jpisdt(nbdysegs) = nbdybeg 1120 jpisft(nbdysegs) = nbdyend 1121 ! 1122 CASE( 'E' ) 1123 IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1 1124 nbdyind = jpiglo - 2 ! set boundary to whole side of model domain. 1125 nbdybeg = 2 1126 nbdyend = jpjglo - 1 1127 ENDIF 1128 nbdysege = nbdysege + 1 1129 npckge(nbdysege) = kb_bdy ! Save bdy package number 1130 jpieob(nbdysege) = nbdyind 1131 jpjedt(nbdysege) = nbdybeg 1132 jpjeft(nbdysege) = nbdyend 1133 ! 1134 CASE( 'W' ) 1135 IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1 1136 nbdyind = 2 ! set boundary to whole side of model domain. 1137 nbdybeg = 2 1138 nbdyend = jpjglo - 1 1139 ENDIF 1140 nbdysegw = nbdysegw + 1 1141 npckgw(nbdysegw) = kb_bdy ! Save bdy package number 1142 jpiwob(nbdysegw) = nbdyind 1143 jpjwdt(nbdysegw) = nbdybeg 1144 jpjwft(nbdysegw) = nbdyend 1145 ! 1146 CASE DEFAULT ; CALL ctl_stop( 'ctypebdy must be N, S, E or W' ) 1147 END SELECT 1148 1149 ! For simplicity we assume that in case of straight bdy, arrays have the same length 1150 ! (even if it is true that last tangential velocity points 1151 ! are useless). This simplifies a little bit boundary data format (and agrees with format 1152 ! used so far in obc package) 1153 1154 knblendta(1:jpbgrd) = (nbdyend - nbdybeg + 1) * nn_rimwidth(kb_bdy) 1155 1156 END SUBROUTINE bdy_read_seg 1157 1158 1378 1159 SUBROUTINE bdy_ctl_seg 1379 1160 !!---------------------------------------------------------------------- … … 1719 1500 END SUBROUTINE bdy_ctl_seg 1720 1501 1721 1502 1503 SUBROUTINE bdy_coords_seg( nbidta, nbjdta, nbrdta ) 1504 !!---------------------------------------------------------------------- 1505 !! *** ROUTINE bdy_coords_seg *** 1506 !! 1507 !! ** Purpose : build nbidta, nbidta, nbrdta for bdy built with segments 1508 !! 1509 !! ** Method : 1510 !! 1511 !!---------------------------------------------------------------------- 1512 INTEGER, DIMENSION(:,:,:), intent( out) :: nbidta, nbjdta, nbrdta ! Index arrays: i and j indices of bdy dta 1513 !! 1514 INTEGER :: ii, ij, ir, iseg 1515 INTEGER :: igrd ! grid type (t=1, u=2, v=3) 1516 INTEGER :: icount ! 1517 INTEGER :: ib_bdy ! bdy number 1518 !!---------------------------------------------------------------------- 1519 1520 ! East 1521 !----- 1522 DO iseg = 1, nbdysege 1523 ib_bdy = npckge(iseg) 1524 ! 1525 ! ------------ T points ------------- 1526 igrd=1 1527 icount=0 1528 DO ir = 1, nn_rimwidth(ib_bdy) 1529 DO ij = jpjedt(iseg), jpjeft(iseg) 1530 icount = icount + 1 1531 nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir 1532 nbjdta(icount, igrd, ib_bdy) = ij 1533 nbrdta(icount, igrd, ib_bdy) = ir 1534 ENDDO 1535 ENDDO 1536 ! 1537 ! ------------ U points ------------- 1538 igrd=2 1539 icount=0 1540 DO ir = 1, nn_rimwidth(ib_bdy) 1541 DO ij = jpjedt(iseg), jpjeft(iseg) 1542 icount = icount + 1 1543 nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 1 - ir 1544 nbjdta(icount, igrd, ib_bdy) = ij 1545 nbrdta(icount, igrd, ib_bdy) = ir 1546 ENDDO 1547 ENDDO 1548 ! 1549 ! ------------ V points ------------- 1550 igrd=3 1551 icount=0 1552 DO ir = 1, nn_rimwidth(ib_bdy) 1553 ! DO ij = jpjedt(iseg), jpjeft(iseg) - 1 1554 DO ij = jpjedt(iseg), jpjeft(iseg) 1555 icount = icount + 1 1556 nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir 1557 nbjdta(icount, igrd, ib_bdy) = ij 1558 nbrdta(icount, igrd, ib_bdy) = ir 1559 ENDDO 1560 nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 1561 nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 1562 ENDDO 1563 ENDDO 1564 ! 1565 ! West 1566 !----- 1567 DO iseg = 1, nbdysegw 1568 ib_bdy = npckgw(iseg) 1569 ! 1570 ! ------------ T points ------------- 1571 igrd=1 1572 icount=0 1573 DO ir = 1, nn_rimwidth(ib_bdy) 1574 DO ij = jpjwdt(iseg), jpjwft(iseg) 1575 icount = icount + 1 1576 nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 1577 nbjdta(icount, igrd, ib_bdy) = ij 1578 nbrdta(icount, igrd, ib_bdy) = ir 1579 ENDDO 1580 ENDDO 1581 ! 1582 ! ------------ U points ------------- 1583 igrd=2 1584 icount=0 1585 DO ir = 1, nn_rimwidth(ib_bdy) 1586 DO ij = jpjwdt(iseg), jpjwft(iseg) 1587 icount = icount + 1 1588 nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 1589 nbjdta(icount, igrd, ib_bdy) = ij 1590 nbrdta(icount, igrd, ib_bdy) = ir 1591 ENDDO 1592 ENDDO 1593 ! 1594 ! ------------ V points ------------- 1595 igrd=3 1596 icount=0 1597 DO ir = 1, nn_rimwidth(ib_bdy) 1598 ! DO ij = jpjwdt(iseg), jpjwft(iseg) - 1 1599 DO ij = jpjwdt(iseg), jpjwft(iseg) 1600 icount = icount + 1 1601 nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 1602 nbjdta(icount, igrd, ib_bdy) = ij 1603 nbrdta(icount, igrd, ib_bdy) = ir 1604 ENDDO 1605 nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 1606 nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 1607 ENDDO 1608 ENDDO 1609 ! 1610 ! North 1611 !----- 1612 DO iseg = 1, nbdysegn 1613 ib_bdy = npckgn(iseg) 1614 ! 1615 ! ------------ T points ------------- 1616 igrd=1 1617 icount=0 1618 DO ir = 1, nn_rimwidth(ib_bdy) 1619 DO ii = jpindt(iseg), jpinft(iseg) 1620 icount = icount + 1 1621 nbidta(icount, igrd, ib_bdy) = ii 1622 nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir 1623 nbrdta(icount, igrd, ib_bdy) = ir 1624 ENDDO 1625 ENDDO 1626 ! 1627 ! ------------ U points ------------- 1628 igrd=2 1629 icount=0 1630 DO ir = 1, nn_rimwidth(ib_bdy) 1631 ! DO ii = jpindt(iseg), jpinft(iseg) - 1 1632 DO ii = jpindt(iseg), jpinft(iseg) 1633 icount = icount + 1 1634 nbidta(icount, igrd, ib_bdy) = ii 1635 nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir 1636 nbrdta(icount, igrd, ib_bdy) = ir 1637 ENDDO 1638 nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 1639 nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 1640 ENDDO 1641 ! 1642 ! ------------ V points ------------- 1643 igrd=3 1644 icount=0 1645 DO ir = 1, nn_rimwidth(ib_bdy) 1646 DO ii = jpindt(iseg), jpinft(iseg) 1647 icount = icount + 1 1648 nbidta(icount, igrd, ib_bdy) = ii 1649 nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 1 - ir 1650 nbrdta(icount, igrd, ib_bdy) = ir 1651 ENDDO 1652 ENDDO 1653 ENDDO 1654 ! 1655 ! South 1656 !----- 1657 DO iseg = 1, nbdysegs 1658 ib_bdy = npckgs(iseg) 1659 ! 1660 ! ------------ T points ------------- 1661 igrd=1 1662 icount=0 1663 DO ir = 1, nn_rimwidth(ib_bdy) 1664 DO ii = jpisdt(iseg), jpisft(iseg) 1665 icount = icount + 1 1666 nbidta(icount, igrd, ib_bdy) = ii 1667 nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 1668 nbrdta(icount, igrd, ib_bdy) = ir 1669 ENDDO 1670 ENDDO 1671 ! 1672 ! ------------ U points ------------- 1673 igrd=2 1674 icount=0 1675 DO ir = 1, nn_rimwidth(ib_bdy) 1676 ! DO ii = jpisdt(iseg), jpisft(iseg) - 1 1677 DO ii = jpisdt(iseg), jpisft(iseg) 1678 icount = icount + 1 1679 nbidta(icount, igrd, ib_bdy) = ii 1680 nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 1681 nbrdta(icount, igrd, ib_bdy) = ir 1682 ENDDO 1683 nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 1684 nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 1685 ENDDO 1686 ! 1687 ! ------------ V points ------------- 1688 igrd=3 1689 icount=0 1690 DO ir = 1, nn_rimwidth(ib_bdy) 1691 DO ii = jpisdt(iseg), jpisft(iseg) 1692 icount = icount + 1 1693 nbidta(icount, igrd, ib_bdy) = ii 1694 nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 1695 nbrdta(icount, igrd, ib_bdy) = ir 1696 ENDDO 1697 ENDDO 1698 ENDDO 1699 1700 1701 END SUBROUTINE bdy_coords_seg 1702 1703 1722 1704 SUBROUTINE bdy_ctl_corn( ib1, ib2 ) 1723 1705 !!---------------------------------------------------------------------- -
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdytides.F90
r11048 r11223 70 70 INTEGER :: inum, igrd 71 71 INTEGER, DIMENSION(3) :: ilen0 !: length of boundary data (from OBC arrays) 72 INTEGER, POINTER, DIMENSION(:) :: nblen, nblenrim ! short cuts73 72 INTEGER :: ios ! Local integer output status for namelist read 74 73 CHARACTER(len=80) :: clfile !: full file name for tidal input file … … 77 76 !! 78 77 TYPE(TIDES_DATA), POINTER :: td !: local short cut 79 TYPE(MAP_POINTER), DIMENSION(jpbgrd) :: ibmap_ptr !: array of pointers to nbmap80 78 !! 81 79 NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta, ln_bdytide_conj 82 80 !!---------------------------------------------------------------------- 83 81 ! 84 IF (nb_bdy>0) THEN 85 IF(lwp) WRITE(numout,*) 86 IF(lwp) WRITE(numout,*) 'bdytide_init : initialization of tidal harmonic forcing at open boundaries' 87 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 88 ENDIF 82 IF(lwp) WRITE(numout,*) 83 IF(lwp) WRITE(numout,*) 'bdytide_init : initialization of tidal harmonic forcing at open boundaries' 84 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 89 85 90 86 REWIND(numnam_cfg) … … 94 90 ! 95 91 td => tides(ib_bdy) 96 nblen => idx_bdy(ib_bdy)%nblen97 nblenrim => idx_bdy(ib_bdy)%nblenrim98 92 99 93 ! Namelist nambdy_tide : tidal harmonic forcing at open boundaries 100 94 filtide(:) = '' 101 95 102 ! Don't REWIND here - may need to read more than one of these namelists.96 REWIND( numnam_ref ) 103 97 READ ( numnam_ref, nambdy_tide, IOSTAT = ios, ERR = 901) 104 98 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy_tide in reference namelist', lwp ) 99 ! Don't REWIND here - may need to read more than one of these namelists. 105 100 READ ( numnam_cfg, nambdy_tide, IOSTAT = ios, ERR = 902 ) 106 101 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nambdy_tide in configuration namelist', lwp ) … … 125 120 ! JC: If FRS scheme is used, we assume that tidal is needed over the whole 126 121 ! relaxation area 127 IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN ; ilen0(:) = nblen (:)128 ELSE ; ilen0(:) = nblenrim(:)122 IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN ; ilen0(:) = idx_bdy(ib_bdy)%nblen (:) 123 ELSE ; ilen0(:) = idx_bdy(ib_bdy)%nblenrim(:) 129 124 ENDIF 130 125 … … 210 205 ALLOCATE( dta_read( MAXVAL(ilen0(1:3)), 1, 1 ) ) 211 206 ! 212 ! Set map structure213 ibmap_ptr(1)%ptr => idx_bdy(ib_bdy)%nbmap(:,1) ; ibmap_ptr(1)%ll_unstruc = ln_coords_file(ib_bdy)214 ibmap_ptr(2)%ptr => idx_bdy(ib_bdy)%nbmap(:,2) ; ibmap_ptr(2)%ll_unstruc = ln_coords_file(ib_bdy)215 ibmap_ptr(3)%ptr => idx_bdy(ib_bdy)%nbmap(:,3) ; ibmap_ptr(3)%ll_unstruc = ln_coords_file(ib_bdy)216 217 207 ! Open files and read in tidal forcing data 218 208 ! ----------------------------------------- … … 222 212 clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_T.nc' 223 213 CALL iom_open( clfile, inum ) 224 CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, ibmap_ptr(1) )214 CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) 225 215 td%ssh0(:,itide,1) = dta_read(1:ilen0(1),1,1) 226 CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1, ibmap_ptr(1) )216 CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) 227 217 td%ssh0(:,itide,2) = dta_read(1:ilen0(1),1,1) 228 218 CALL iom_close( inum ) … … 230 220 clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_U.nc' 231 221 CALL iom_open( clfile, inum ) 232 CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, i bmap_ptr(2) )222 CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) 233 223 td%u0(:,itide,1) = dta_read(1:ilen0(2),1,1) 234 CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, i bmap_ptr(2) )224 CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) 235 225 td%u0(:,itide,2) = dta_read(1:ilen0(2),1,1) 236 226 CALL iom_close( inum ) … … 238 228 clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_V.nc' 239 229 CALL iom_open( clfile, inum ) 240 CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, i bmap_ptr(3) )230 CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) 241 231 td%v0(:,itide,1) = dta_read(1:ilen0(3),1,1) 242 CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, i bmap_ptr(3) )232 CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) 243 233 td%v0(:,itide,2) = dta_read(1:ilen0(3),1,1) 244 234 CALL iom_close( inum ) … … 272 262 273 263 274 SUBROUTINE bdytide_update( kt, idx, dta, td, jit, time_offset )264 SUBROUTINE bdytide_update( kt, idx, dta, td, kit, kt_offset ) 275 265 !!---------------------------------------------------------------------- 276 266 !! *** SUBROUTINE bdytide_update *** … … 283 273 TYPE(OBC_DATA) , INTENT(inout) :: dta ! OBC external data 284 274 TYPE(TIDES_DATA) , INTENT(inout) :: td ! tidal harmonics data 285 INTEGER, OPTIONAL, INTENT(in ) :: jit ! Barotropic timestep counter (for timesplitting option)286 INTEGER, OPTIONAL, INTENT(in ) :: time_offset ! time offset in units of timesteps. NB. if jit275 INTEGER, OPTIONAL, INTENT(in ) :: kit ! Barotropic timestep counter (for timesplitting option) 276 INTEGER, OPTIONAL, INTENT(in ) :: kt_offset ! time offset in units of timesteps. NB. if kit 287 277 ! ! is present then units = subcycle timesteps. 288 ! ! time_offset = 0 => get data at "now" time level289 ! ! time_offset = -1 => get data at "before" time level290 ! ! time_offset = +1 => get data at "after" time level278 ! ! kt_offset = 0 => get data at "now" time level 279 ! ! kt_offset = -1 => get data at "before" time level 280 ! ! kt_offset = +1 => get data at "after" time level 291 281 ! ! etc. 292 282 ! … … 303 293 304 294 zflag=1 305 IF ( PRESENT( jit) ) THEN306 IF ( jit /= 1 ) zflag=0295 IF ( PRESENT(kit) ) THEN 296 IF ( kit /= 1 ) zflag=0 307 297 ENDIF 308 298 … … 323 313 324 314 time_add = 0 325 IF( PRESENT( time_offset) ) THEN326 time_add = time_offset315 IF( PRESENT(kt_offset) ) THEN 316 time_add = kt_offset 327 317 ENDIF 328 318 329 IF( PRESENT( jit) ) THEN330 z_arg = ((kt-kt_tide) * rdt + ( jit+0.5_wp*(time_add-1)) * rdt / REAL(nn_baro,wp) )319 IF( PRESENT(kit) ) THEN 320 z_arg = ((kt-kt_tide) * rdt + (kit+0.5_wp*(time_add-1)) * rdt / REAL(nn_baro,wp) ) 331 321 ELSE 332 322 z_arg = ((kt-kt_tide)+time_add) * rdt … … 361 351 362 352 363 SUBROUTINE bdy_dta_tides( kt, kit, time_offset )353 SUBROUTINE bdy_dta_tides( kt, kit, kt_offset ) 364 354 !!---------------------------------------------------------------------- 365 355 !! *** SUBROUTINE bdy_dta_tides *** … … 370 360 INTEGER, INTENT(in) :: kt ! Main timestep counter 371 361 INTEGER, OPTIONAL, INTENT(in) :: kit ! Barotropic timestep counter (for timesplitting option) 372 INTEGER, OPTIONAL, INTENT(in) :: time_offset! time offset in units of timesteps. NB. if kit362 INTEGER, OPTIONAL, INTENT(in) :: kt_offset ! time offset in units of timesteps. NB. if kit 373 363 ! ! is present then units = subcycle timesteps. 374 ! ! time_offset = 0 => get data at "now" time level375 ! ! time_offset = -1 => get data at "before" time level376 ! ! time_offset = +1 => get data at "after" time level364 ! ! kt_offset = 0 => get data at "now" time level 365 ! ! kt_offset = -1 => get data at "before" time level 366 ! ! kt_offset = +1 => get data at "after" time level 377 367 ! ! etc. 378 368 ! … … 389 379 390 380 time_add = 0 391 IF( PRESENT( time_offset) ) THEN392 time_add = time_offset381 IF( PRESENT(kt_offset) ) THEN 382 time_add = kt_offset 393 383 ENDIF 394 384 … … 435 425 ! If time splitting, initialize arrays from slow varying open boundary data: 436 426 IF ( PRESENT(kit) ) THEN 437 IF ( dta_bdy(ib_bdy)%l l_ssh) dta_bdy(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy_s(ib_bdy)%ssh(1:ilen0(1))438 IF ( dta_bdy(ib_bdy)%l l_u2d ) dta_bdy(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy_s(ib_bdy)%u2d(1:ilen0(2))439 IF ( dta_bdy(ib_bdy)%l l_v2d ) dta_bdy(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy_s(ib_bdy)%v2d(1:ilen0(3))427 IF ( dta_bdy(ib_bdy)%lneed_ssh ) dta_bdy(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy_s(ib_bdy)%ssh(1:ilen0(1)) 428 IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) dta_bdy(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy_s(ib_bdy)%u2d(1:ilen0(2)) 429 IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) dta_bdy(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy_s(ib_bdy)%v2d(1:ilen0(3)) 440 430 ENDIF 441 431 ! … … 447 437 z_sist = zramp * SIN( z_sarg ) 448 438 ! 449 IF ( dta_bdy(ib_bdy)%l l_ssh ) THEN439 IF ( dta_bdy(ib_bdy)%lneed_ssh ) THEN 450 440 igrd=1 ! SSH on tracer grid 451 441 DO ib = 1, ilen0(igrd) … … 456 446 ENDIF 457 447 ! 458 IF ( dta_bdy(ib_bdy)%l l_u2d ) THEN448 IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) THEN 459 449 igrd=2 ! U grid 460 450 DO ib = 1, ilen0(igrd) … … 463 453 & tides(ib_bdy)%u(ib,itide,2)*z_sist ) 464 454 END DO 465 ENDIF466 !467 IF ( dta_bdy(ib_bdy)%ll_v2d ) THEN468 455 igrd=3 ! V grid 469 456 DO ib = 1, ilen0(igrd) -
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/DOM/dommsk.F90
r10425 r11223 101 101 & cn_ice, nn_ice_dta, & 102 102 & rn_ice_tem, rn_ice_sal, rn_ice_age, & 103 & ln_vol, nn_volctl, nn_rimwidth , nb_jpk_bdy103 & ln_vol, nn_volctl, nn_rimwidth 104 104 !!--------------------------------------------------------------------- 105 105 ! -
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/DYN/dynspg_ts.F90
r10742 r11223 718 718 ! ! ------------------ 719 719 ! Update only tidal forcing at open boundaries 720 IF( ln_bdy .AND. ln_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset= noffset+1 )721 IF( ln_tide_pot .AND. ln_tide ) CALL upd_tide ( kt, kit=jn, time_offset= noffset )720 IF( ln_bdy .AND. ln_tide ) CALL bdy_dta_tides( kt, kit=jn, kt_offset= noffset+1 ) 721 IF( ln_tide_pot .AND. ln_tide ) CALL upd_tide ( kt, kit=jn, kt_offset= noffset ) 722 722 ! 723 723 ! Set extrapolation coefficients for predictor step: -
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/IOM/iom.F90
r11194 r11223 835 835 836 836 837 FUNCTION iom_varid ( kiomid, cdvar, kdimsz, kndims, ld stop )837 FUNCTION iom_varid ( kiomid, cdvar, kdimsz, kndims, lduld, ldstop ) 838 838 !!----------------------------------------------------------------------- 839 839 !! *** FUNCTION iom_varid *** … … 844 844 CHARACTER(len=*) , INTENT(in ) :: cdvar ! name of the variable 845 845 INTEGER, DIMENSION(:), INTENT( out), OPTIONAL :: kdimsz ! size of each dimension 846 INTEGER, INTENT( out), OPTIONAL :: kndims ! size of the dimensions 846 INTEGER , INTENT( out), OPTIONAL :: kndims ! number of dimensions 847 LOGICAL , INTENT( out), OPTIONAL :: lduld ! true if the last dimension is unlimited (time) 847 848 LOGICAL , INTENT(in ), OPTIONAL :: ldstop ! stop if looking for non-existing variable (default = .TRUE.) 848 849 ! … … 874 875 iiv = iiv + 1 875 876 IF( iiv <= jpmax_vars ) THEN 876 iom_varid = iom_nf90_varid( kiomid, cdvar, iiv, kdimsz, kndims )877 iom_varid = iom_nf90_varid( kiomid, cdvar, iiv, kdimsz, kndims, lduld ) 877 878 ELSE 878 879 CALL ctl_stop( trim(clinfo), 'Too many variables in the file '//iom_file(kiomid)%name, & … … 892 893 ENDIF 893 894 IF( PRESENT(kndims) ) kndims = iom_file(kiomid)%ndims(iiv) 895 IF( PRESENT( lduld) ) lduld = iom_file(kiomid)%luld( iiv) 894 896 ENDIF 895 897 ENDIF -
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/IOM/iom_nf90.F90
r10522 r11223 187 187 188 188 189 FUNCTION iom_nf90_varid ( kiomid, cdvar, kiv, kdimsz, kndims )189 FUNCTION iom_nf90_varid ( kiomid, cdvar, kiv, kdimsz, kndims, lduld ) 190 190 !!----------------------------------------------------------------------- 191 191 !! *** FUNCTION iom_varid *** … … 198 198 INTEGER, DIMENSION(:), INTENT( out), OPTIONAL :: kdimsz ! size of the dimensions 199 199 INTEGER, INTENT( out), OPTIONAL :: kndims ! size of the dimensions 200 LOGICAL , INTENT( out), OPTIONAL :: lduld ! true if the last dimension is unlimited (time) 200 201 ! 201 202 INTEGER :: iom_nf90_varid ! iom variable Id … … 251 252 ENDIF 252 253 IF( PRESENT(kndims) ) kndims = iom_file(kiomid)%ndims(kiv) 254 IF( PRESENT( lduld) ) lduld = iom_file(kiomid)%luld(kiv) 253 255 ELSE 254 256 iom_nf90_varid = -1 ! variable not found, return error code: -1 -
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/LBC/mppini.F90
r11192 r11223 167 167 & cn_ice, nn_ice_dta, & 168 168 & rn_ice_tem, rn_ice_sal, rn_ice_age, & 169 & ln_vol, nn_volctl, nn_rimwidth , nb_jpk_bdy169 & ln_vol, nn_volctl, nn_rimwidth 170 170 !!---------------------------------------------------------------------- 171 171 -
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/SBC/fldread.F90
r10425 r11223 80 80 INTEGER :: nreclast ! last record to be read in the current file 81 81 CHARACTER(len = 256) :: lsmname ! current name of the NetCDF mask file acting as a key 82 ! ! variables related to BDY 82 83 INTEGER :: igrd ! grid type for bdy data 83 84 INTEGER :: ibdy ! bdy set id number 85 INTEGER, POINTER, DIMENSION(:) :: imap ! Array of integer pointers to 1D arrays 86 LOGICAL :: ltotvel ! total velocity or not (T/F) 84 87 END TYPE FLD 85 86 TYPE, PUBLIC :: MAP_POINTER !: Map from input data file to local domain87 INTEGER, POINTER, DIMENSION(:) :: ptr ! Array of integer pointers to 1D arrays88 LOGICAL :: ll_unstruc ! Unstructured (T) or structured (F) boundary data file89 END TYPE MAP_POINTER90 88 91 89 !$AGRIF_DO_NOT_TREAT … … 129 127 CONTAINS 130 128 131 SUBROUTINE fld_read( kt, kn_fsbc, sd, map, kit, kt_offset, jpk_bdy, fvl)129 SUBROUTINE fld_read( kt, kn_fsbc, sd, kit, kt_offset ) 132 130 !!--------------------------------------------------------------------- 133 131 !! *** ROUTINE fld_read *** … … 144 142 INTEGER , INTENT(in ) :: kn_fsbc ! sbc computation period (in time step) 145 143 TYPE(FLD), INTENT(inout), DIMENSION(:) :: sd ! input field related variables 146 TYPE(MAP_POINTER),INTENT(in), OPTIONAL, DIMENSION(:) :: map ! global-to-local mapping indices147 144 INTEGER , INTENT(in ), OPTIONAL :: kit ! subcycle timestep for timesplitting option 148 145 INTEGER , INTENT(in ), OPTIONAL :: kt_offset ! provide fields at time other than "now" … … 150 147 ! ! kt_offset = +1 => fields at "after" time level 151 148 ! ! etc. 152 INTEGER , INTENT(in ), OPTIONAL :: jpk_bdy ! number of vertical levels in the BDY data153 LOGICAL , INTENT(in ), OPTIONAL :: fvl ! number of vertical levels in the BDY data154 149 !! 155 150 INTEGER :: itmp ! local variable … … 166 161 REAL(wp) :: ztintb ! ratio applied to before records when doing time interpolation 167 162 CHARACTER(LEN=1000) :: clfmt ! write format 168 TYPE(MAP_POINTER) :: imap ! global-to-local mapping indices169 163 !!--------------------------------------------------------------------- 170 164 ll_firstcall = kt == nit000 … … 175 169 ENDIF 176 170 IF( PRESENT(kt_offset) ) it_offset = kt_offset 177 178 imap%ptr => NULL()179 171 180 172 ! Note that shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar … … 188 180 IF( ll_firstcall ) THEN ! initialization 189 181 DO jf = 1, imf 190 IF( PRESENT(map) ) imap = map(jf) 191 IF( PRESENT(jpk_bdy) ) THEN 192 CALL fld_init( kn_fsbc, sd(jf), imap, jpk_bdy, fvl ) ! read each before field (put them in after as they will be swapped) 193 ELSE 194 CALL fld_init( kn_fsbc, sd(jf), imap ) ! read each before field (put them in after as they will be swapped) 195 ENDIF 182 IF( TRIM(sd(jf)%clrootname) == 'NOT USED' ) CYCLE 183 CALL fld_init( kn_fsbc, sd(jf) ) ! read each before field (put them in after as they will be swapped) 196 184 END DO 197 185 IF( lwp ) CALL wgt_print() ! control print … … 202 190 ! 203 191 DO jf = 1, imf ! --- loop over field --- ! 204 192 193 IF( TRIM(sd(jf)%clrootname) == 'NOT USED' ) CYCLE 194 205 195 IF( isecsbc > sd(jf)%nrec_a(2) .OR. ll_firstcall ) THEN ! read/update the after data? 206 207 IF( PRESENT(map) ) imap = map(jf) ! temporary definition of map208 196 209 197 sd(jf)%nrec_b(:) = sd(jf)%nrec_a(:) ! swap before record informations … … 222 210 itmp = sd(jf)%nrec_a(1) ! temporary storage 223 211 sd(jf)%nrec_a(1) = sd(jf)%nreclast ! read the last record of the file currently opened 224 CALL fld_get( sd(jf) , imap )! read after data212 CALL fld_get( sd(jf) ) ! read after data 225 213 sd(jf)%fdta(:,:,:,1) = sd(jf)%fdta(:,:,:,2) ! re-swap before record field 226 214 sd(jf)%nrec_b(1) = sd(jf)%nrec_a(1) ! update before record informations … … 240 228 & .AND. sd(jf)%nrec_b(1) /= sd(jf)%nrec_a(1) - 1 ) THEN 241 229 sd(jf)%nrec_a(1) = sd(jf)%nrec_a(1) - 1 ! move back to before record 242 CALL fld_get( sd(jf) , imap )! read after data230 CALL fld_get( sd(jf) ) ! read after data 243 231 sd(jf)%fdta(:,:,:,1) = sd(jf)%fdta(:,:,:,2) ! re-swap before record field 244 232 sd(jf)%nrec_b(1) = sd(jf)%nrec_a(1) ! update before record informations … … 285 273 286 274 ! read after data 287 IF( PRESENT(jpk_bdy) ) THEN 288 CALL fld_get( sd(jf), imap, jpk_bdy, fvl) 289 ELSE 290 CALL fld_get( sd(jf), imap ) 291 ENDIF 275 CALL fld_get( sd(jf) ) 276 292 277 ENDIF ! read new data? 293 278 END DO ! --- end loop over field --- ! … … 296 281 297 282 DO jf = 1, imf ! --- loop over field --- ! 283 ! 284 IF( TRIM(sd(jf)%clrootname) == 'NOT USED' ) CYCLE 298 285 ! 299 286 IF( sd(jf)%ln_tint ) THEN ! temporal interpolation … … 327 314 328 315 329 SUBROUTINE fld_init( kn_fsbc, sdjf , map , jpk_bdy, fvl)316 SUBROUTINE fld_init( kn_fsbc, sdjf ) 330 317 !!--------------------------------------------------------------------- 331 318 !! *** ROUTINE fld_init *** … … 336 323 INTEGER , INTENT(in ) :: kn_fsbc ! sbc computation period (in time step) 337 324 TYPE(FLD), INTENT(inout) :: sdjf ! input field related variables 338 TYPE(MAP_POINTER),INTENT(in) :: map ! global-to-local mapping indices339 INTEGER , INTENT(in), OPTIONAL :: jpk_bdy ! number of vertical levels in the BDY data340 LOGICAL , INTENT(in), OPTIONAL :: fvl ! number of vertical levels in the BDY data341 325 !! 342 326 LOGICAL :: llprevyr ! are we reading previous year file? … … 351 335 CHARACTER(LEN=1000) :: clfmt ! write format 352 336 !!--------------------------------------------------------------------- 337 ! 353 338 llprevyr = .FALSE. 354 339 llprevmth = .FALSE. … … 433 418 ! 434 419 ! read before data in after arrays(as we will swap it later) 435 IF( PRESENT(jpk_bdy) ) THEN 436 CALL fld_get( sdjf, map, jpk_bdy, fvl ) 437 ELSE 438 CALL fld_get( sdjf, map ) 439 ENDIF 420 CALL fld_get( sdjf ) 440 421 ! 441 422 clfmt = "(' fld_init : time-interpolation for ', a, ' read previous record = ', i6, ' at time = ', f7.2, ' days')" … … 613 594 614 595 615 SUBROUTINE fld_get( sdjf , map, jpk_bdy, fvl)596 SUBROUTINE fld_get( sdjf ) 616 597 !!--------------------------------------------------------------------- 617 598 !! *** ROUTINE fld_get *** … … 620 601 !!---------------------------------------------------------------------- 621 602 TYPE(FLD) , INTENT(inout) :: sdjf ! input field related variables 622 TYPE(MAP_POINTER), INTENT(in ) :: map ! global-to-local mapping indices623 INTEGER , INTENT(in), OPTIONAL :: jpk_bdy ! number of vertical levels in the bdy data624 LOGICAL , INTENT(in), OPTIONAL :: fvl ! number of vertical levels in the bdy data625 603 ! 626 604 INTEGER :: ipk ! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk ) … … 634 612 ipk = SIZE( sdjf%fnow, 3 ) 635 613 ! 636 IF( ASSOCIATED(map%ptr) ) THEN 637 IF( PRESENT(jpk_bdy) ) THEN 638 IF( sdjf%ln_tint ) THEN ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2), & 639 sdjf%nrec_a(1), map, sdjf%igrd, sdjf%ibdy, jpk_bdy, fvl ) 640 ELSE ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,: ), & 641 sdjf%nrec_a(1), map, sdjf%igrd, sdjf%ibdy, jpk_bdy, fvl ) 642 ENDIF 643 ELSE 644 IF( sdjf%ln_tint ) THEN ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1), map ) 645 ELSE ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,: ), sdjf%nrec_a(1), map ) 646 ENDIF 647 ENDIF 614 IF( ASSOCIATED(sdjf%imap) ) THEN 615 IF( sdjf%ln_tint ) THEN ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2), & 616 & sdjf%nrec_a(1), sdjf%imap, sdjf%igrd, sdjf%ibdy, sdjf%ltotvel ) 617 ELSE ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,: ), & 618 & sdjf%nrec_a(1), sdjf%imap, sdjf%igrd, sdjf%ibdy, sdjf%ltotvel ) 619 ENDIF 648 620 ELSE IF( LEN(TRIM(sdjf%wgtname)) > 0 ) THEN 649 621 CALL wgt_list( sdjf, iw ) … … 700 672 END SUBROUTINE fld_get 701 673 702 SUBROUTINE fld_map( num, clvar, dta, nrec, map, igrd, ibdy, jpk_bdy, fvl ) 674 675 SUBROUTINE fld_map( knum, cdvar, pdta, krec, kmap, kgrd, kbdy, ldtotvel ) 703 676 !!--------------------------------------------------------------------- 704 677 !! *** ROUTINE fld_map *** … … 707 680 !! using a general mapping (for open boundaries) 708 681 !!---------------------------------------------------------------------- 709 710 USE bdy_oce, ONLY: ln_bdy, idx_bdy, dta_global, dta_global_z, dta_global_dz, dta_global2, dta_global2_z, dta_global2_dz ! workspace to read in global data arrays 711 712 INTEGER , INTENT(in ) :: num ! stream number 713 CHARACTER(LEN=*) , INTENT(in ) :: clvar ! variable name 714 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: dta ! output field on model grid (2 dimensional) 715 INTEGER , INTENT(in ) :: nrec ! record number to read (ie time slice) 716 TYPE(MAP_POINTER) , INTENT(in ) :: map ! global-to-local mapping indices 717 INTEGER , INTENT(in), OPTIONAL :: igrd, ibdy, jpk_bdy ! grid type, set number and number of vertical levels in the bdy data 718 LOGICAL , INTENT(in), OPTIONAL :: fvl ! grid type, set number and number of vertical levels in the bdy data 719 INTEGER :: jpkm1_bdy! number of vertical levels in the bdy data minus 1 720 !! 721 INTEGER :: ipi ! length of boundary data on local process 722 INTEGER :: ipj ! length of dummy dimension ( = 1 ) 723 INTEGER :: ipk ! number of vertical levels of dta ( 2D: ipk=1 ; 3D: ipk=jpk ) 724 INTEGER :: ilendta ! length of data in file 725 INTEGER :: idvar ! variable ID 726 INTEGER :: ib, ik, ji, jj ! loop counters 727 INTEGER :: ierr 728 REAL(wp) :: fv ! fillvalue 729 REAL(wp), POINTER, DIMENSION(:,:,:) :: dta_read ! work space for global data 730 REAL(wp), POINTER, DIMENSION(:,:,:) :: dta_read_z ! work space for global data 731 REAL(wp), POINTER, DIMENSION(:,:,:) :: dta_read_dz ! work space for global data 732 !!--------------------------------------------------------------------- 733 ! 734 ipi = SIZE( dta, 1 ) 735 ipj = 1 736 ipk = SIZE( dta, 3 ) 737 ! 738 idvar = iom_varid( num, clvar ) 739 ilendta = iom_file(num)%dimsz(1,idvar) 740 741 IF ( ln_bdy ) THEN 742 ipj = iom_file(num)%dimsz(2,idvar) 743 IF( map%ll_unstruc) THEN ! unstructured open boundary data file 744 dta_read => dta_global 745 IF( PRESENT(jpk_bdy) ) THEN 746 IF( jpk_bdy>0 ) THEN 747 dta_read_z => dta_global_z 748 dta_read_dz => dta_global_dz 749 jpkm1_bdy = jpk_bdy-1 750 ENDIF 751 ENDIF 752 ELSE ! structured open boundary file 753 dta_read => dta_global2 754 IF( PRESENT(jpk_bdy) ) THEN 755 IF( jpk_bdy>0 ) THEN 756 dta_read_z => dta_global2_z 757 dta_read_dz => dta_global2_dz 758 jpkm1_bdy = jpk_bdy-1 759 ENDIF 760 ENDIF 761 ENDIF 762 ENDIF 763 764 IF(lwp) WRITE(numout,*) 'Dim size for ', TRIM(clvar),' is ', ilendta 765 IF(lwp) WRITE(numout,*) 'Number of levels for ',TRIM(clvar),' is ', ipk 766 ! 767 SELECT CASE( ipk ) 768 CASE(1) ; 769 CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1 ), nrec ) 770 IF ( map%ll_unstruc) THEN ! unstructured open boundary data file 771 DO ib = 1, ipi 772 DO ik = 1, ipk 773 dta(ib,1,ik) = dta_read(map%ptr(ib),1,ik) 774 END DO 775 END DO 776 ELSE ! we assume that this is a structured open boundary file 777 DO ib = 1, ipi 778 jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 779 ji=map%ptr(ib)-(jj-1)*ilendta 780 DO ik = 1, ipk 781 dta(ib,1,ik) = dta_read(ji,jj,ik) 782 END DO 783 END DO 784 ENDIF 682 INTEGER , INTENT(in ) :: knum ! stream number 683 CHARACTER(LEN=*) , INTENT(in ) :: cdvar ! variable name 684 REAL(wp), DIMENSION(:,:,:), INTENT( out) :: pdta ! bdy output field on model grid 685 INTEGER , INTENT(in ) :: krec ! record number to read (ie time slice) 686 INTEGER , DIMENSION(:) , INTENT(in ) :: kmap ! global-to-local bdy mapping indices 687 ! optional variables used for vertical interpolation: 688 INTEGER, OPTIONAL , INTENT(in ) :: kgrd ! grid type (t, u, v) 689 INTEGER, OPTIONAL , INTENT(in ) :: kbdy ! bdy number 690 LOGICAL, OPTIONAL , INTENT(in ) :: ldtotvel ! true if toal ( = barotrop + barocline) velocity 691 !! 692 INTEGER :: ipi ! length of boundary data on local process 693 INTEGER :: ipj ! length of dummy dimension ( = 1 ) 694 INTEGER :: ipk ! number of vertical levels of pdta ( 2D: ipk=1 ; 3D: ipk=jpk ) 695 INTEGER :: ipkb ! number of vertical levels in boundary data file 696 INTEGER :: idvar ! variable ID 697 INTEGER :: indims ! number of dimensions of the variable 698 INTEGER, DIMENSION(4) :: idimsz ! size of variable dimensions 699 REAL(wp) :: zfv ! fillvalue 700 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zz_read ! work space for global boundary data 701 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zdta_read ! work space local data requiring vertical interpolation 702 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zdta_read_z ! work space local data requiring vertical interpolation 703 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zdta_read_dz ! work space local data requiring vertical interpolation 704 CHARACTER(LEN=1),DIMENSION(3) :: clgrid 705 LOGICAL :: lluld ! is the variable using the unlimited dimension 706 !!--------------------------------------------------------------------- 707 ! 708 clgrid = (/'t','u','v'/) 709 ! 710 ipi = SIZE( pdta, 1 ) 711 ipj = SIZE( pdta, 2 ) ! must be equal to 1 712 ipk = SIZE( pdta, 3 ) 713 ! 714 idvar = iom_varid( knum, cdvar, kndims = indims, kdimsz = idimsz, lduld = lluld ) 715 IF( indims == 4 .OR. ( indims == 3 .AND. .NOT. lluld ) ) THEN ; ipkb = idimsz(3) ! xy(zl)t or xy(zl) 716 ELSE ; ipkb = 1 ! xy or xyt 717 ENDIF 718 ! 719 ALLOCATE( zz_read( idimsz(1), idimsz(2), ipkb ) ) ! ++++++++ !!! this can be very big... 720 ! 721 IF( ipk == 1 ) THEN 722 723 IF( ipkb /= 1 ) CALL ctl_stop( 'fld_map : we must have ipkb = 1 to read surface data' ) 724 CALL iom_get ( knum, jpdom_unknown, cdvar, zz_read(:,:,1), krec ) ! call iom_get with a 2D file 725 CALL fld_map_core( zz_read, kmap, pdta ) 785 726 786 727 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 787 728 ! Do we include something here to adjust barotropic velocities ! 788 729 ! in case of a depth difference between bdy files and ! 789 ! bathymetry in the case ln_ full_vel = .false. and jpk_bdy>0?!730 ! bathymetry in the case ln_totvel = .false. and ipkb>0? ! 790 731 ! [as the enveloping and parital cells could change H] ! 791 732 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 792 733 793 CASE DEFAULT ; 794 795 IF( PRESENT(jpk_bdy) .AND. jpk_bdy>0 ) THEN ! boundary data not on model grid: vertical interpolation 796 CALL iom_getatt(num, '_FillValue', fv, cdvar=clvar ) 797 dta_read(:,:,:) = -ABS(fv) 798 dta_read_z(:,:,:) = 0._wp 799 dta_read_dz(:,:,:) = 0._wp 800 CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:jpk_bdy), nrec ) 801 SELECT CASE( igrd ) 802 CASE(1) 803 CALL iom_get ( num, jpdom_unknown, 'gdept', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) ) 804 CALL iom_get ( num, jpdom_unknown, 'e3t', dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) ) 805 CASE(2) 806 CALL iom_get ( num, jpdom_unknown, 'gdepu', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) ) 807 CALL iom_get ( num, jpdom_unknown, 'e3u', dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) ) 808 CASE(3) 809 CALL iom_get ( num, jpdom_unknown, 'gdepv', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) ) 810 CALL iom_get ( num, jpdom_unknown, 'e3v', dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) ) 811 END SELECT 812 813 IF ( ln_bdy ) & 814 CALL fld_bdy_interp(dta_read, dta_read_z, dta_read_dz, map, jpk_bdy, igrd, ibdy, fv, dta, fvl, ilendta) 815 816 ELSE ! boundary data assumed to be on model grid 817 818 CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:ipk), nrec ) 819 IF ( map%ll_unstruc) THEN ! unstructured open boundary data file 820 DO ib = 1, ipi 821 DO ik = 1, ipk 822 dta(ib,1,ik) = dta_read(map%ptr(ib),1,ik) 823 END DO 734 ELSE 735 ! 736 CALL iom_get ( knum, jpdom_unknown, cdvar, zz_read(:,:,:), krec ) ! call iom_get with a 3D file 737 ! 738 IF( ipkb /= ipk ) THEN ! boundary data not on model vertical grid : vertical interpolation 739 ! 740 IF( ipk == jpk .AND. iom_varid(knum,'gdep'//clgrid(kgrd)) /= -1 .AND. iom_varid(knum,'e3'//clgrid(kgrd)) /= -1 ) THEN 741 742 ALLOCATE( zdta_read(ipi,ipj,ipkb), zdta_read_z(ipi,ipj,ipkb), zdta_read_dz(ipi,ipj,ipkb) ) 743 744 CALL fld_map_core( zz_read, kmap, zdta_read ) 745 CALL iom_get ( knum, jpdom_unknown, 'gdep'//clgrid(kgrd), zz_read ) ! read only once? Potential temporal evolution? 746 CALL fld_map_core( zz_read, kmap, zdta_read_z ) 747 CALL iom_get ( knum, jpdom_unknown, 'e3'//clgrid(kgrd), zz_read ) ! read only once? Potential temporal evolution? 748 CALL fld_map_core( zz_read, kmap, zdta_read_dz ) 749 750 CALL iom_getatt(knum, '_FillValue', zfv, cdvar=cdvar ) 751 CALL fld_bdy_interp(zdta_read, zdta_read_z, zdta_read_dz, pdta, kgrd, kbdy, zfv, ldtotvel) 752 DEALLOCATE( zdta_read, zdta_read_z, zdta_read_dz ) 753 754 ELSE 755 IF( ipk /= jpk ) CALL ctl_stop( 'fld_map : this should be an impossible case...' ) 756 WRITE(ctmp1,*) 'fld_map : vertical interpolation for bdy variable '//TRIM(cdvar)//' requires ' 757 IF( iom_varid(knum, 'gdep'//clgrid(kgrd)) == -1 ) CALL ctl_stop( ctmp1//'gdep'//clgrid(kgrd)//' variable' ) 758 IF( iom_varid(knum, 'e3'//clgrid(kgrd)) == -1 ) CALL ctl_stop( ctmp1// 'e3'//clgrid(kgrd)//' variable' ) 759 760 ENDIF 761 ! 762 ELSE ! bdy data assumed to be the same levels as bdy variables 763 ! 764 CALL fld_map_core( zz_read, kmap, pdta ) 765 ! 766 ENDIF ! ipkb /= ipk 767 ENDIF ! ipk == 1 768 769 DEALLOCATE( zz_read ) 770 771 END SUBROUTINE fld_map 772 773 774 SUBROUTINE fld_map_core( pdta_read, kmap, pdta_bdy ) 775 !!--------------------------------------------------------------------- 776 !! *** ROUTINE fld_map_core *** 777 !! 778 !! ** Purpose : inner core of fld_map 779 !!---------------------------------------------------------------------- 780 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pdta_read ! global boundary data 781 INTEGER, DIMENSION(: ), INTENT(in ) :: kmap ! global-to-local bdy mapping indices 782 REAL(wp), DIMENSION(:,:,:), INTENT( out) :: pdta_bdy ! bdy output field on model grid 783 !! 784 INTEGER, DIMENSION(3) :: idim_read, idim_bdy ! arrays dimensions 785 INTEGER :: ji, jj, jk, jb ! loop counters 786 INTEGER :: im1 787 !!--------------------------------------------------------------------- 788 ! 789 idim_read = SHAPE( pdta_read ) 790 idim_bdy = SHAPE( pdta_bdy ) 791 ! 792 ! in all cases: idim_bdy(2) == 1 .AND. idim_read(1) * idim_read(2) == idim_bdy(1) 793 ! structured BDY with rimwidth > 1 : idim_read(2) == rimwidth /= 1 794 ! structured BDY with rimwidth == 1 or unstructured BDY: idim_read(2) == 1 795 ! 796 IF( idim_read(2) > 1 ) THEN ! structured BDY with rimwidth > 1 797 DO jk = 1, idim_bdy(3) 798 DO jb = 1, idim_bdy(1) 799 im1 = kmap(jb) - 1 800 jj = im1 / idim_read(1) + 1 801 ji = MOD( im1, idim_read(1) ) + 1 802 pdta_bdy(jb,1,jk) = pdta_read(ji,jj,jk) 824 803 END DO 825 ELSE ! we assume that this is a structured open boundary file 826 DO ib = 1, ipi 827 jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 828 ji=map%ptr(ib)-(jj-1)*ilendta 829 DO ik = 1, ipk 830 dta(ib,1,ik) = dta_read(ji,jj,ik) 831 END DO 804 END DO 805 ELSE 806 DO jk = 1, idim_bdy(3) 807 DO jb = 1, idim_bdy(1) ! horizontal remap of bdy data on the local bdy 808 pdta_bdy(jb,1,jk) = pdta_read(kmap(jb),1,jk) 832 809 END DO 833 ENDIF 834 ENDIF ! PRESENT(jpk_bdy) 835 END SELECT 836 837 END SUBROUTINE fld_map 810 END DO 811 ENDIF 812 813 END SUBROUTINE fld_map_core 838 814 839 SUBROUTINE fld_bdy_interp(dta_read, dta_read_z, dta_read_dz, map, jpk_bdy, igrd, ibdy, fv, dta, fvl, ilendta)840 815 816 SUBROUTINE fld_bdy_interp(pdta_read, pdta_read_z, pdta_read_dz, pdta, kgrd, kbdy, pfv, ldtotvel) 841 817 !!--------------------------------------------------------------------- 842 818 !! *** ROUTINE fld_bdy_interp *** … … 847 823 USE bdy_oce, ONLY: idx_bdy ! indexing for map <-> ij transformation 848 824 849 REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in ) :: dta_read ! work space for global data 850 REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in ) :: dta_read_z ! work space for global data 851 REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in ) :: dta_read_dz ! work space for global data 852 REAL(wp) , INTENT(in) :: fv ! fillvalue and alternative -ABS(fv) 853 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: dta ! output field on model grid (2 dimensional) 854 TYPE(MAP_POINTER) , INTENT(in ) :: map ! global-to-local mapping indices 855 LOGICAL , INTENT(in), OPTIONAL :: fvl ! grid type, set number and number of vertical levels in the bdy data 856 INTEGER , INTENT(in) :: igrd, ibdy, jpk_bdy ! number of levels in bdy data 857 INTEGER , INTENT(in) :: ilendta ! length of data in file 858 !! 859 INTEGER :: ipi ! length of boundary data on local process 860 INTEGER :: ipj ! length of dummy dimension ( = 1 ) 861 INTEGER :: ipk ! number of vertical levels of dta ( 2D: ipk=1 ; 3D: ipk=jpk ) 862 INTEGER :: jpkm1_bdy ! number of levels in bdy data minus 1 863 INTEGER :: ib, ik, ikk ! loop counters 864 INTEGER :: ji, jj, zij, zjj ! temporary indices 865 REAL(wp) :: zl, zi, zh ! tmp variable for current depth and interpolation factor 866 REAL(wp) :: fv_alt, ztrans, ztrans_new ! fillvalue and alternative -ABS(fv) 867 CHARACTER (LEN=10) :: ibstr 825 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pdta_read ! data read in bdy file 826 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pdta_read_z ! depth of the data read in bdy file 827 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pdta_read_dz ! thickness of the levels in bdy file 828 REAL(wp), DIMENSION(:,:,:), INTENT( out) :: pdta ! output field on model grid (2 dimensional) 829 REAL(wp) , INTENT(in ) :: pfv ! fillvalue of the data read in bdy file 830 LOGICAL , INTENT(in ) :: ldtotvel ! true if toal ( = barotrop + barocline) velocity 831 INTEGER , INTENT(in ) :: kgrd ! grid type (t, u, v) 832 INTEGER , INTENT(in ) :: kbdy ! bdy number 833 !! 834 INTEGER :: ipi ! length of boundary data on local process 835 INTEGER :: ipkb ! number of vertical levels in boundary data file 836 INTEGER :: jb, ji, jj, jk, jkb ! loop counters 837 REAL(wp) :: zcoef 838 REAL(wp) :: zl, zi, zh ! tmp variable for current depth and interpolation factor 839 REAL(wp) :: zfv_alt, ztrans, ztrans_new ! fillvalue and alternative -ABS(pfv) 840 REAL(wp), DIMENSION(jpk) :: zdepth, zdhalf ! level and half-level depth 868 841 !!--------------------------------------------------------------------- 869 842 870 871 ipi = SIZE( dta, 1 ) 872 ipj = SIZE( dta_read, 2 ) 873 ipk = SIZE( dta, 3 ) 874 jpkm1_bdy = jpk_bdy-1 843 ipi = SIZE( pdta, 1 ) 844 ipkb = SIZE( pdta_read, 3 ) 875 845 876 fv_alt = -ABS(fv) ! set _FillValue < 0 as we make use of MAXVAL and MAXLOC later 877 DO ib = 1, ipi 878 zij = idx_bdy(ibdy)%nbi(ib,igrd) 879 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 880 IF(narea==2) WRITE(*,*) 'MAPI', ib, igrd, map%ptr(ib), narea-1, zij, zjj 881 ENDDO 882 ! 883 IF ( map%ll_unstruc ) THEN ! unstructured open boundary data file 884 885 DO ib = 1, ipi 886 DO ik = 1, jpk_bdy 887 IF( ( dta_read(map%ptr(ib),1,ik) == fv ) ) THEN 888 dta_read_z(map%ptr(ib),1,ik) = fv_alt ! safety: put fillvalue into external depth field so consistent with data 889 dta_read_dz(map%ptr(ib),1,ik) = 0._wp ! safety: put 0._wp into external thickness factors to ensure transport is correct 846 zfv_alt = -ABS(pfv) ! set _FillValue < 0 as we make use of MAXVAL and MAXLOC later 847 ! 848 WHERE( pdta_read == pfv ) 849 pdta_read_z = zfv_alt ! safety: put fillvalue into external depth field so consistent with data 850 pdta_read_dz = 0._wp ! safety: put 0._wp into external thickness factors to ensure transport is correct 851 ENDWHERE 852 853 DO jb = 1, ipi 854 ji = idx_bdy(kbdy)%nbi(jb,kgrd) 855 jj = idx_bdy(kbdy)%nbj(jb,kgrd) 856 zh = SUM(pdta_read_dz(jb,1,:) ) 857 ! 858 ! Warnings to flag differences in the input and model topgraphy - is this useful/necessary? 859 SELECT CASE( kgrd ) 860 CASE(1) 861 IF( ABS( (zh - ht_n(ji,jj)) / ht_n(ji,jj)) * tmask(ji,jj,1) > 0.01_wp ) THEN 862 WRITE(ctmp1,"(I10.10)") jb 863 CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ctmp1)//' by more than 1%') 864 ! IF(lwp) WRITE(numout,*) 'DEPTHT', zh, sum(e3t_n(ji,jj,:), mask=tmask(ji,jj,:)==1), ht_n(ji,jj), jb, jb, ji, jj 865 ENDIF 866 CASE(2) 867 IF( ABS( (zh - hu_n(ji,jj)) * r1_hu_n(ji,jj)) * umask(ji,jj,1) > 0.01_wp ) THEN 868 WRITE(ctmp1,"(I10.10)") jb 869 CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ctmp1)//' by more than 1%') 870 ! IF(lwp) WRITE(numout,*) 'DEPTHU', zh, SUM(e3u_n(ji,jj,:), mask=umask(ji,jj,:)==1), SUM(umask(ji,jj,:)), & 871 ! & hu_n(ji,jj), jb, jb, ji, jj, narea-1, pdta_read(jb,1,:) 872 ENDIF 873 CASE(3) 874 IF( ABS( (zh - hv_n(ji,jj)) * r1_hv_n(ji,jj)) * vmask(ji,jj,1) > 0.01_wp ) THEN 875 WRITE(ctmp1,"(I10.10)") jb 876 CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ctmp1)//' by more than 1%') 877 ENDIF 878 END SELECT 879 ! 880 SELECT CASE( kgrd ) 881 CASE(1) 882 ! depth of T points: 883 zdepth(:) = gdept_n(ji,jj,:) 884 CASE(2) 885 ! depth of U points: we must not use gdept_n as we don't want to do a communication 886 ! --> copy what is done for gdept_n in domvvl... 887 zdhalf(1) = 0.0_wp 888 zdepth(1) = 0.5_wp * e3uw_n(ji,jj,1) 889 DO jk = 2, jpk ! vertical sum 890 ! zcoef = umask - wumask ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 891 ! ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 892 ! ! 0.5 where jk = mikt 893 !!gm ??????? BUG ? gdept_n as well as gde3w_n does not include the thickness of ISF ?? 894 zcoef = ( umask(ji,jj,jk) - wumask(ji,jj,jk) ) 895 zdhalf(jk) = zdhalf(jk-1) + e3u_n(ji,jj,jk-1) 896 zdepth(jk) = zcoef * ( zdhalf(jk ) + 0.5 * e3uw_n(ji,jj,jk)) & 897 & + (1-zcoef) * ( zdepth(jk-1) + e3uw_n(ji,jj,jk)) 898 END DO 899 CASE(3) 900 ! depth of V points: we must not use gdept_n as we don't want to do a communication 901 ! --> copy what is done for gdept_n in domvvl... 902 zdhalf(1) = 0.0_wp 903 zdepth(1) = 0.5_wp * e3vw_n(ji,jj,1) 904 DO jk = 2, jpk ! vertical sum 905 ! zcoef = vmask - wvmask ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 906 ! ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 907 ! ! 0.5 where jk = mikt 908 !!gm ??????? BUG ? gdept_n as well as gde3w_n does not include the thickness of ISF ?? 909 zcoef = ( vmask(ji,jj,jk) - wvmask(ji,jj,jk) ) 910 zdhalf(jk) = zdhalf(jk-1) + e3v_n(ji,jj,jk-1) 911 zdepth(jk) = zcoef * ( zdhalf(jk ) + 0.5 * e3vw_n(ji,jj,jk)) & 912 & + (1-zcoef) * ( zdepth(jk-1) + e3vw_n(ji,jj,jk)) 913 END DO 914 END SELECT 915 ! 916 DO jk = 1, jpk 917 IF( zdepth(jk) < pdta_read_z(jb,1, 1) ) THEN ! above the first level of external data 918 pdta(jb,1,jk) = pdta_read(jb,1,1) 919 ELSEIF( zdepth(jk) > pdta_read_z(jb,1,ipkb) ) THEN ! below the last level of external data 920 pdta(jb,1,jk) = pdta_read(jb,1,MAXLOC(pdta_read_z(jb,1,:),1)) 921 ELSE ! inbetween: vertical interpolation between jkb & jkb+1 922 DO jkb = 1, ipkb-1 ! when gdept_n(jkb) < zdepth(jk) < gdept_n(jkb+1) 923 IF( ( ( zdepth(jk) - pdta_read_z(jb,1,jkb) ) * ( zdepth(jk) - pdta_read_z(jb,1,jkb+1) ) <= 0._wp ) & 924 & .AND. ( pdta_read_z(jb,1,jkb+1) /= zfv_alt) ) THEN ! linear interpolation between 2 levels 925 zi = ( zdepth(jk) - pdta_read_z(jb,1,jkb) ) / ( pdta_read_z(jb,1,jkb+1) - pdta_read_z(jb,1,jkb) ) 926 pdta(jb,1,jk) = pdta_read(jb,1,jkb) + ( pdta_read (jb,1,jkb+1) - pdta_read (jb,1,jkb) ) * zi 927 ENDIF 928 END DO 929 ENDIF 930 END DO ! jpk 931 ! 932 END DO ! ipi 933 934 IF(kgrd == 2) THEN ! do we need to adjust the transport term? 935 DO jb = 1, ipi 936 ji = idx_bdy(kbdy)%nbi(jb,kgrd) 937 jj = idx_bdy(kbdy)%nbj(jb,kgrd) 938 zh = SUM(pdta_read_dz(jb,1,:) ) 939 ztrans = 0._wp 940 ztrans_new = 0._wp 941 DO jkb = 1, ipkb ! calculate transport on input grid 942 ztrans = ztrans + pdta_read(jb,1,jkb) * pdta_read_dz(jb, 1,jkb) 943 ENDDO 944 DO jk = 1, jpk ! calculate transport on model grid 945 ztrans_new = ztrans_new + pdta(jb,1,jk ) * e3u_n(ji,jj,jk ) * umask(ji,jj,jk) 946 ENDDO 947 DO jk = 1, jpk ! make transport correction 948 IF(ldtotvel) THEN ! bdy data are total velocity so adjust bt transport term to match input data 949 pdta(jb,1,jk) = ( pdta(jb,1,jk) + ( ztrans - ztrans_new ) * r1_hu_n(ji,jj) ) * umask(ji,jj,jk) 950 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 951 pdta(jb,1,jk) = pdta(jb,1,jk) + ( 0._wp - ztrans_new ) * r1_hu_n(ji,jj) * umask(ji,jj,jk) 890 952 ENDIF 891 953 ENDDO 892 ENDDO 893 894 DO ib = 1, ipi 895 zij = idx_bdy(ibdy)%nbi(ib,igrd) 896 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 897 zh = SUM(dta_read_dz(map%ptr(ib),1,:) ) 898 ! Warnings to flag differences in the input and model topgraphy - is this useful/necessary? 899 SELECT CASE( igrd ) 900 CASE(1) 901 IF( ABS( (zh - ht_n(zij,zjj)) / ht_n(zij,zjj)) * tmask(zij,zjj,1) > 0.01_wp ) THEN 902 WRITE(ibstr,"(I10.10)") map%ptr(ib) 903 CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 904 ! IF(lwp) WRITE(*,*) 'DEPTHT', zh, sum(e3t_n(zij,zjj,:), mask=tmask(zij,zjj,:)==1), ht_n(zij,zjj), map%ptr(ib), ib, zij, zjj 905 ENDIF 906 CASE(2) 907 IF( ABS( (zh - hu_n(zij,zjj)) * r1_hu_n(zij,zjj)) * umask(zij,zjj,1) > 0.01_wp ) THEN 908 WRITE(ibstr,"(I10.10)") map%ptr(ib) 909 CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 910 IF(lwp) WRITE(*,*) 'DEPTHU', zh, sum(e3u_n(zij,zjj,:), mask=umask(zij,zjj,:)==1), sum(umask(zij,zjj,:)), & 911 & hu_n(zij,zjj), map%ptr(ib), ib, zij, zjj, narea-1 , & 912 & dta_read(map%ptr(ib),1,:) 913 ENDIF 914 CASE(3) 915 IF( ABS( (zh - hv_n(zij,zjj)) * r1_hv_n(zij,zjj)) * vmask(zij,zjj,1) > 0.01_wp ) THEN 916 WRITE(ibstr,"(I10.10)") map%ptr(ib) 917 CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 918 ENDIF 919 END SELECT 920 DO ik = 1, ipk 921 SELECT CASE( igrd ) 922 CASE(1) 923 zl = gdept_n(zij,zjj,ik) ! if using in step could use fsdept instead of gdept_n? 924 CASE(2) 925 IF(ln_sco) THEN 926 zl = ( gdept_n(zij,zjj,ik) + gdept_n(zij+1,zjj,ik) ) * 0.5_wp ! if using in step could use fsdept instead of gdept_n? 927 ELSE 928 zl = MIN( gdept_n(zij,zjj,ik), gdept_n(zij+1,zjj,ik) ) 929 ENDIF 930 CASE(3) 931 IF(ln_sco) THEN 932 zl = ( gdept_n(zij,zjj,ik) + gdept_n(zij,zjj+1,ik) ) * 0.5_wp ! if using in step could use fsdept instead of gdept_n? 933 ELSE 934 zl = MIN( gdept_n(zij,zjj,ik), gdept_n(zij,zjj+1,ik) ) 935 ENDIF 936 END SELECT 937 IF( zl < dta_read_z(map%ptr(ib),1,1) ) THEN ! above the first level of external data 938 dta(ib,1,ik) = dta_read(map%ptr(ib),1,1) 939 ELSEIF( zl > MAXVAL(dta_read_z(map%ptr(ib),1,:),1) ) THEN ! below the last level of external data 940 dta(ib,1,ik) = dta_read(map%ptr(ib),1,MAXLOC(dta_read_z(map%ptr(ib),1,:),1)) 941 ELSE ! inbetween : vertical interpolation between ikk & ikk+1 942 DO ikk = 1, jpkm1_bdy ! when gdept_n(ikk) < zl < gdept_n(ikk+1) 943 IF( ( (zl-dta_read_z(map%ptr(ib),1,ikk)) * (zl-dta_read_z(map%ptr(ib),1,ikk+1)) <= 0._wp) & 944 & .AND. (dta_read_z(map%ptr(ib),1,ikk+1) /= fv_alt)) THEN 945 zi = ( zl - dta_read_z(map%ptr(ib),1,ikk) ) / & 946 & ( dta_read_z(map%ptr(ib),1,ikk+1) - dta_read_z(map%ptr(ib),1,ikk) ) 947 dta(ib,1,ik) = dta_read(map%ptr(ib),1,ikk) + & 948 & ( dta_read(map%ptr(ib),1,ikk+1) - dta_read(map%ptr(ib),1,ikk) ) * zi 949 ENDIF 950 END DO 951 ENDIF 952 END DO 953 END DO 954 955 IF(igrd == 2) THEN ! do we need to adjust the transport term? 956 DO ib = 1, ipi 957 zij = idx_bdy(ibdy)%nbi(ib,igrd) 958 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 959 zh = SUM(dta_read_dz(map%ptr(ib),1,:) ) 960 ztrans = 0._wp 961 ztrans_new = 0._wp 962 DO ik = 1, jpk_bdy ! calculate transport on input grid 963 ztrans = ztrans + dta_read(map%ptr(ib),1,ik) * dta_read_dz(map%ptr(ib),1,ik) 964 ENDDO 965 DO ik = 1, ipk ! calculate transport on model grid 966 ztrans_new = ztrans_new + dta(ib,1,ik) * e3u_n(zij,zjj,ik) * umask(zij,zjj,ik) 967 ENDDO 968 DO ik = 1, ipk ! make transport correction 969 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 970 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hu_n(zij,zjj) ) * umask(zij,zjj,ik) 971 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 972 IF( ABS(ztrans * r1_hu_n(zij,zjj)) > 0.01_wp ) & 973 & CALL ctl_warn('fld_bdy_interp: barotropic component of > 0.01 ms-1 found in baroclinic velocities at') 974 dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hu_n(zij,zjj) * umask(zij,zjj,ik) 975 ENDIF 976 ENDDO 954 ENDDO 955 ENDIF 956 957 IF(kgrd == 3) THEN ! do we need to adjust the transport term? 958 DO jb = 1, ipi 959 ji = idx_bdy(kbdy)%nbi(jb,kgrd) 960 jj = idx_bdy(kbdy)%nbj(jb,kgrd) 961 zh = SUM(pdta_read_dz(jb,1,:) ) 962 ztrans = 0._wp 963 ztrans_new = 0._wp 964 DO jkb = 1, ipkb ! calculate transport on input grid 965 ztrans = ztrans + pdta_read(jb,1,jkb) * pdta_read_dz(jb, 1,jkb) 977 966 ENDDO 978 ENDIF 979 980 IF(igrd == 3) THEN ! do we need to adjust the transport term? 981 DO ib = 1, ipi 982 zij = idx_bdy(ibdy)%nbi(ib,igrd) 983 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 984 zh = SUM(dta_read_dz(map%ptr(ib),1,:) ) 985 ztrans = 0._wp 986 ztrans_new = 0._wp 987 DO ik = 1, jpk_bdy ! calculate transport on input grid 988 ztrans = ztrans + dta_read(map%ptr(ib),1,ik) * dta_read_dz(map%ptr(ib),1,ik) 989 ENDDO 990 DO ik = 1, ipk ! calculate transport on model grid 991 ztrans_new = ztrans_new + dta(ib,1,ik) * e3v_n(zij,zjj,ik) * vmask(zij,zjj,ik) 992 ENDDO 993 DO ik = 1, ipk ! make transport correction 994 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 995 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hv_n(zij,zjj) ) * vmask(zij,zjj,ik) 996 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 997 dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hv_n(zij,zjj) * vmask(zij,zjj,ik) 998 ENDIF 999 ENDDO 967 DO jk = 1, jpk ! calculate transport on model grid 968 ztrans_new = ztrans_new + pdta(jb,1,jk ) * e3v_n(ji,jj,jk ) * vmask(ji,jj,jk) 1000 969 ENDDO 1001 ENDIF 1002 1003 ELSE ! structured open boundary file 1004 1005 DO ib = 1, ipi 1006 jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 1007 ji=map%ptr(ib)-(jj-1)*ilendta 1008 DO ik = 1, jpk_bdy 1009 IF( ( dta_read(ji,jj,ik) == fv ) ) THEN 1010 dta_read_z(ji,jj,ik) = fv_alt ! safety: put fillvalue into external depth field so consistent with data 1011 dta_read_dz(ji,jj,ik) = 0._wp ! safety: put 0._wp into external thickness factors to ensure transport is correct 970 DO jk = 1, jpk ! make transport correction 971 IF(ldtotvel) THEN ! bdy data are total velocity so adjust bt transport term to match input data 972 pdta(jb,1,jk) = ( pdta(jb,1,jk) + ( ztrans - ztrans_new ) * r1_hv_n(ji,jj) ) * vmask(ji,jj,jk) 973 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 974 pdta(jb,1,jk) = pdta(jb,1,jk) + ( 0._wp - ztrans_new ) * r1_hv_n(ji,jj) * vmask(ji,jj,jk) 1012 975 ENDIF 1013 976 ENDDO 1014 ENDDO 1015 1016 1017 DO ib = 1, ipi 1018 jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 1019 ji=map%ptr(ib)-(jj-1)*ilendta 1020 zij = idx_bdy(ibdy)%nbi(ib,igrd) 1021 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 1022 zh = SUM(dta_read_dz(ji,jj,:) ) 1023 ! Warnings to flag differences in the input and model topgraphy - is this useful/necessary? 1024 SELECT CASE( igrd ) 1025 CASE(1) 1026 IF( ABS( (zh - ht_n(zij,zjj)) / ht_n(zij,zjj)) * tmask(zij,zjj,1) > 0.01_wp ) THEN 1027 WRITE(ibstr,"(I10.10)") map%ptr(ib) 1028 CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 1029 ! IF(lwp) WRITE(*,*) 'DEPTHT', zh, sum(e3t_n(zij,zjj,:), mask=tmask(zij,zjj,:)==1), ht_n(zij,zjj), map%ptr(ib), ib, zij, zjj 1030 ENDIF 1031 CASE(2) 1032 IF( ABS( (zh - hu_n(zij,zjj)) * r1_hu_n(zij,zjj)) * umask(zij,zjj,1) > 0.01_wp ) THEN 1033 WRITE(ibstr,"(I10.10)") map%ptr(ib) 1034 CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 1035 ENDIF 1036 CASE(3) 1037 IF( ABS( (zh - hv_n(zij,zjj)) * r1_hv_n(zij,zjj)) * vmask(zij,zjj,1) > 0.01_wp ) THEN 1038 WRITE(ibstr,"(I10.10)") map%ptr(ib) 1039 CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 1040 ENDIF 1041 END SELECT 1042 DO ik = 1, ipk 1043 SELECT CASE( igrd ) ! coded for sco - need zco and zps option using min 1044 CASE(1) 1045 zl = gdept_n(zij,zjj,ik) ! if using in step could use fsdept instead of gdept_n? 1046 CASE(2) 1047 IF(ln_sco) THEN 1048 zl = ( gdept_n(zij,zjj,ik) + gdept_n(zij+1,zjj,ik) ) * 0.5_wp ! if using in step could use fsdept instead of gdept_n? 1049 ELSE 1050 zl = MIN( gdept_n(zij,zjj,ik), gdept_n(zij+1,zjj,ik) ) 1051 ENDIF 1052 CASE(3) 1053 IF(ln_sco) THEN 1054 zl = ( gdept_n(zij,zjj,ik) + gdept_n(zij,zjj+1,ik) ) * 0.5_wp ! if using in step could use fsdept instead of gdept_n? 1055 ELSE 1056 zl = MIN( gdept_n(zij,zjj,ik), gdept_n(zij,zjj+1,ik) ) 1057 ENDIF 1058 END SELECT 1059 IF( zl < dta_read_z(ji,jj,1) ) THEN ! above the first level of external data 1060 dta(ib,1,ik) = dta_read(ji,jj,1) 1061 ELSEIF( zl > MAXVAL(dta_read_z(ji,jj,:),1) ) THEN ! below the last level of external data 1062 dta(ib,1,ik) = dta_read(ji,jj,MAXLOC(dta_read_z(ji,jj,:),1)) 1063 ELSE ! inbetween : vertical interpolation between ikk & ikk+1 1064 DO ikk = 1, jpkm1_bdy ! when gdept_n(ikk) < zl < gdept_n(ikk+1) 1065 IF( ( (zl-dta_read_z(ji,jj,ikk)) * (zl-dta_read_z(ji,jj,ikk+1)) <= 0._wp) & 1066 & .AND. (dta_read_z(ji,jj,ikk+1) /= fv_alt)) THEN 1067 zi = ( zl - dta_read_z(ji,jj,ikk) ) / & 1068 & ( dta_read_z(ji,jj,ikk+1) - dta_read_z(ji,jj,ikk) ) 1069 dta(ib,1,ik) = dta_read(ji,jj,ikk) + & 1070 & ( dta_read(ji,jj,ikk+1) - dta_read(ji,jj,ikk) ) * zi 1071 ENDIF 1072 END DO 1073 ENDIF 1074 END DO 1075 END DO 1076 1077 IF(igrd == 2) THEN ! do we need to adjust the transport term? 1078 DO ib = 1, ipi 1079 jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 1080 ji=map%ptr(ib)-(jj-1)*ilendta 1081 zij = idx_bdy(ibdy)%nbi(ib,igrd) 1082 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 1083 zh = SUM(dta_read_dz(ji,jj,:) ) 1084 ztrans = 0._wp 1085 ztrans_new = 0._wp 1086 DO ik = 1, jpk_bdy ! calculate transport on input grid 1087 ztrans = ztrans + dta_read(ji,jj,ik) * dta_read_dz(ji,jj,ik) 1088 ENDDO 1089 DO ik = 1, ipk ! calculate transport on model grid 1090 ztrans_new = ztrans_new + dta(ib,1,ik) * e3u_n(zij,zjj,ik) * umask(zij,zjj,ik) 1091 ENDDO 1092 DO ik = 1, ipk ! make transport correction 1093 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 1094 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hu_n(zij,zjj) ) * umask(zij,zjj,ik) 1095 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 1096 dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hu_n(zij,zjj) ) * umask(zij,zjj,ik) 1097 ENDIF 1098 ENDDO 1099 ENDDO 1100 ENDIF 1101 1102 IF(igrd == 3) THEN ! do we need to adjust the transport term? 1103 DO ib = 1, ipi 1104 jj = 1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 1105 ji = map%ptr(ib)-(jj-1)*ilendta 1106 zij = idx_bdy(ibdy)%nbi(ib,igrd) 1107 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 1108 zh = SUM(dta_read_dz(ji,jj,:) ) 1109 ztrans = 0._wp 1110 ztrans_new = 0._wp 1111 DO ik = 1, jpk_bdy ! calculate transport on input grid 1112 ztrans = ztrans + dta_read(ji,jj,ik) * dta_read_dz(ji,jj,ik) 1113 ENDDO 1114 DO ik = 1, ipk ! calculate transport on model grid 1115 ztrans_new = ztrans_new + dta(ib,1,ik) * e3v_n(zij,zjj,ik) * vmask(zij,zjj,ik) 1116 ENDDO 1117 DO ik = 1, ipk ! make transport correction 1118 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 1119 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hv_n(zij,zjj) ) * vmask(zij,zjj,ik) 1120 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 1121 dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hv_n(zij,zjj) ) * vmask(zij,zjj,ik) 1122 ENDIF 1123 ENDDO 1124 ENDDO 1125 ENDIF 1126 1127 ENDIF ! endif unstructured or structured 1128 977 ENDDO 978 ENDIF 979 1129 980 END SUBROUTINE fld_bdy_interp 1130 981 … … 1151 1002 imf = SIZE( sd ) 1152 1003 DO ju = 1, imf 1004 IF( TRIM(sd(ju)%clrootname) == 'NOT USED' ) CYCLE 1153 1005 ill = LEN_TRIM( sd(ju)%vcomp ) 1154 1006 DO jn = 2-COUNT((/sd(ju)%ln_tint/)), 2 … … 1159 1011 iv = -1 1160 1012 DO jv = 1, imf 1013 IF( TRIM(sd(jv)%clrootname) == 'NOT USED' ) CYCLE 1161 1014 IF( TRIM(sd(jv)%vcomp) == TRIM(clcomp) ) iv = jv 1162 1015 END DO … … 1299 1152 ! 1300 1153 DO jf = 1, SIZE(sdf) 1301 sdf(jf)%clrootname = TRIM( cdir )//TRIM( sdf_n(jf)%clname ) 1154 sdf(jf)%clrootname = sdf_n(jf)%clname 1155 IF( TRIM(sdf_n(jf)%clname) /= 'NOT USED' ) sdf(jf)%clrootname = TRIM( cdir )//sdf(jf)%clrootname 1302 1156 sdf(jf)%clname = "not yet defined" 1303 1157 sdf(jf)%nfreqh = sdf_n(jf)%nfreqh … … 1308 1162 sdf(jf)%num = -1 1309 1163 sdf(jf)%wgtname = " " 1310 IF( LEN( TRIM(sdf_n(jf)%wname) ) > 0 ) sdf(jf)%wgtname = TRIM( cdir )// TRIM( sdf_n(jf)%wname )1164 IF( LEN( TRIM(sdf_n(jf)%wname) ) > 0 ) sdf(jf)%wgtname = TRIM( cdir )//sdf_n(jf)%wname 1311 1165 sdf(jf)%lsmname = " " 1312 IF( LEN( TRIM(sdf_n(jf)%lname) ) > 0 ) sdf(jf)%lsmname = TRIM( cdir )// TRIM( sdf_n(jf)%lname )1166 IF( LEN( TRIM(sdf_n(jf)%lname) ) > 0 ) sdf(jf)%lsmname = TRIM( cdir )//sdf_n(jf)%lname 1313 1167 sdf(jf)%vcomp = sdf_n(jf)%vcomp 1314 1168 sdf(jf)%rotn(:) = .TRUE. ! pretend to be rotated -> won't try to rotate data before the first call to fld_get … … 1317 1171 IF( sdf(jf)%cltype(1:4) == 'week' .AND. sdf(jf)%ln_clim ) & 1318 1172 & CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdf(jf)%clrootname)//') needs ln_clim = .FALSE.') 1319 sdf(jf)%nreclast = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during fld_clopn 1173 sdf(jf)%nreclast = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during fld_clopn 1174 sdf(jf)%igrd = 0 1175 sdf(jf)%ibdy = 0 1176 sdf(jf)%imap => NULL() 1177 sdf(jf)%ltotvel = .FALSE. 1320 1178 END DO 1321 1179 ! -
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/SBC/updtide.F90
r10068 r11223 27 27 CONTAINS 28 28 29 SUBROUTINE upd_tide( kt, kit, time_offset )29 SUBROUTINE upd_tide( kt, kit, kt_offset ) 30 30 !!---------------------------------------------------------------------- 31 31 !! *** ROUTINE upd_tide *** … … 39 39 INTEGER, INTENT(in) :: kt ! ocean time-step index 40 40 INTEGER, INTENT(in), OPTIONAL :: kit ! external mode sub-time-step index (lk_dynspg_ts=T) 41 INTEGER, INTENT(in), OPTIONAL :: time_offset ! time offset in number41 INTEGER, INTENT(in), OPTIONAL :: kt_offset ! time offset in number 42 42 ! of internal steps (lk_dynspg_ts=F) 43 43 ! of external steps (lk_dynspg_ts=T) 44 44 ! 45 INTEGER :: joffset ! local integer45 INTEGER :: ioffset ! local integer 46 46 INTEGER :: ji, jj, jk ! dummy loop indices 47 47 REAL(wp) :: zt, zramp ! local scalar … … 52 52 zt = ( kt - kt_tide ) * rdt 53 53 ! 54 joffset = 055 IF( PRESENT( time_offset ) ) joffset = time_offset54 ioffset = 0 55 IF( PRESENT( kt_offset ) ) ioffset = kt_offset 56 56 ! 57 57 IF( PRESENT( kit ) ) THEN 58 zt = zt + ( kit + joffset - 1 ) * rdt / REAL( nn_baro, wp )58 zt = zt + ( kit + ioffset - 1 ) * rdt / REAL( nn_baro, wp ) 59 59 ELSE 60 zt = zt + joffset * rdt60 zt = zt + ioffset * rdt 61 61 ENDIF 62 62 ! … … 70 70 IF( ln_tide_ramp ) THEN ! linear increase if asked 71 71 zt = ( kt - nit000 ) * rdt 72 IF( PRESENT( kit ) ) zt = zt + ( kit + joffset -1) * rdt / REAL( nn_baro, wp )72 IF( PRESENT( kit ) ) zt = zt + ( kit + ioffset -1) * rdt / REAL( nn_baro, wp ) 73 73 zramp = MIN( MAX( zt / (rdttideramp*rday) , 0._wp ) , 1._wp ) 74 74 pot_astro(:,:) = zramp * pot_astro(:,:) -
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/step.F90
r10364 r11223 112 112 IF( ln_tide ) CALL sbc_tide( kstp ) ! update tide potential 113 113 IF( ln_apr_dyn ) CALL sbc_apr ( kstp ) ! atmospheric pressure (NB: call before bdy_dta which needs ssh_ib) 114 IF( ln_bdy ) CALL bdy_dta ( kstp, time_offset=+1 ) ! update dynamic & tracer data at open boundaries114 IF( ln_bdy ) CALL bdy_dta ( kstp, kt_offset = +1 ) ! update dynamic & tracer data at open boundaries 115 115 CALL sbc ( kstp ) ! Sea Boundary Condition (including sea-ice) 116 116
Note: See TracChangeset
for help on using the changeset viewer.