Changeset 12367
- Timestamp:
- 2020-02-11T19:37:34+01:00 (3 years ago)
- File:
-
- 1 edited
-
NEMO/trunk/src/OCE/SBC/fldread.F90 (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/OCE/SBC/fldread.F90
r12248 r12367 833 833 834 834 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pdta_read ! data read in bdy file 835 REAL(wp), DIMENSION(:,:,:), INTENT(in out) :: pdta_read_z ! depth of the data read in bdy file836 REAL(wp), DIMENSION(:,:,:), INTENT(in out) :: pdta_read_dz ! thickness of the levels in bdy file835 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pdta_read_z ! depth of the data read in bdy file 836 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pdta_read_dz ! thickness of the levels in bdy file 837 837 REAL(wp), DIMENSION(:,:,:), INTENT( out) :: pdta ! output field on model grid (2 dimensional) 838 838 REAL(wp) , INTENT(in ) :: pfv ! fillvalue of the data read in bdy file … … 841 841 INTEGER , INTENT(in ) :: kbdy ! bdy number 842 842 !! 843 INTEGER :: ipi! length of boundary data on local process844 INTEGER :: ipkb! number of vertical levels in boundary data file845 INTEGER :: jb, ji, jj, jk, jkb ! loop counters846 REAL(wp) :: zcoef847 REAL(wp) :: zl, zi, zh ! tmp variable for current depth and interpolation factor848 REAL(wp) :: zfv_alt, ztrans, ztrans_new ! fillvalue and alternative -ABS(pfv)849 REAL(wp), DIMENSION(jpk) :: zdepth, zdhalf! level and half-level depth843 INTEGER :: ipi ! length of boundary data on local process 844 INTEGER :: ipkb ! number of vertical levels in boundary data file 845 INTEGER :: ipkmax ! number of vertical levels in boundary data file where no mask 846 INTEGER :: jb, ji, jj, jk, jkb ! loop counters 847 REAL(wp) :: zcoef, zi ! 848 REAL(wp) :: ztrans, ztrans_new ! transports 849 REAL(wp), DIMENSION(jpk) :: zdepth, zdhalf ! level and half-level depth 850 850 !!--------------------------------------------------------------------- 851 851 … … 853 853 ipkb = SIZE( pdta_read, 3 ) 854 854 855 zfv_alt = -ABS(pfv) ! set _FillValue < 0 as we make use of MAXVAL and MAXLOC later856 !857 WHERE( pdta_read == pfv )858 pdta_read_z = zfv_alt ! safety: put fillvalue into external depth field so consistent with data859 pdta_read_dz = 0._wp ! safety: put 0._wp into external thickness factors to ensure transport is correct860 ENDWHERE861 862 855 DO jb = 1, ipi 863 856 ji = idx_bdy(kbdy)%nbi(jb,kgrd) 864 857 jj = idx_bdy(kbdy)%nbj(jb,kgrd) 865 zh = SUM(pdta_read_dz(jb,1,:) ) 866 ! 858 ! 859 ! --- max jk where input data /= FillValue --- ! 860 ipkmax = 1 861 DO jkb = 2, ipkb 862 IF( pdta_read(jb,1,jkb) /= pfv ) ipkmax = MAX( ipkmax, jkb ) 863 END DO 864 ! 865 ! --- calculate depth at t,u,v points --- ! 867 866 SELECT CASE( kgrd ) 868 CASE(1) 869 ! depth of T points: 867 CASE(1) ! depth of T points: 870 868 zdepth(:) = gdept_n(ji,jj,:) 871 CASE(2) 872 ! depth of U points: we must not use gdept_n as we don't want to do a communication 873 ! --> copy what is done for gdept_n in domvvl... 869 CASE(2) ! depth of U points: we must not use gdept_n as we don't want to do a communication 870 ! --> copy what is done for gdept_n in domvvl... 874 871 zdhalf(1) = 0.0_wp 875 872 zdepth(1) = 0.5_wp * e3uw_n(ji,jj,1) … … 881 878 zcoef = ( umask(ji,jj,jk) - wumask(ji,jj,jk) ) 882 879 zdhalf(jk) = zdhalf(jk-1) + e3u_n(ji,jj,jk-1) 883 zdepth(jk) = zcoef * ( zdhalf(jk ) + 0.5 * e3uw_n(ji,jj,jk)) &884 & + (1 -zcoef) * ( zdepth(jk-1) + e3uw_n(ji,jj,jk))880 zdepth(jk) = zcoef * ( zdhalf(jk ) + 0.5 * e3uw_n(ji,jj,jk)) & 881 & + (1.-zcoef) * ( zdepth(jk-1) + e3uw_n(ji,jj,jk)) 885 882 END DO 886 CASE(3) 887 ! depth of V points: we must not use gdept_n as we don't want to do a communication 888 ! --> copy what is done for gdept_n in domvvl... 883 CASE(3) ! depth of V points: we must not use gdept_n as we don't want to do a communication 884 ! --> copy what is done for gdept_n in domvvl... 889 885 zdhalf(1) = 0.0_wp 890 886 zdepth(1) = 0.5_wp * e3vw_n(ji,jj,1) … … 896 892 zcoef = ( vmask(ji,jj,jk) - wvmask(ji,jj,jk) ) 897 893 zdhalf(jk) = zdhalf(jk-1) + e3v_n(ji,jj,jk-1) 898 zdepth(jk) = zcoef * ( zdhalf(jk ) + 0.5 * e3vw_n(ji,jj,jk)) &899 & + (1 -zcoef) * ( zdepth(jk-1) + e3vw_n(ji,jj,jk))894 zdepth(jk) = zcoef * ( zdhalf(jk ) + 0.5 * e3vw_n(ji,jj,jk)) & 895 & + (1.-zcoef) * ( zdepth(jk-1) + e3vw_n(ji,jj,jk)) 900 896 END DO 901 897 END SELECT 902 ! 903 DO jk = 1, jpk904 IF( zdepth(jk) < pdta_read_z(jb,1, 1) ) THEN ! above the first level of external data905 pdta(jb,1,jk) = pdta_read(jb,1,1)906 ELSEIF( zdepth(jk) > pdta_read_z(jb,1,ipkb) ) THEN ! below the last level of external data907 pdta(jb,1,jk) = pdta_read(jb,1,MAXLOC(pdta_read_z(jb,1,:),1))908 ELSE ! inbetween: vertical interpolation between jkb & jkb+1909 DO jkb = 1, ipkb-1 ! when gdept_n(jkb) < zdepth(jk) < gdept_n(jkb+1)910 IF( ( ( zdepth(jk) - pdta_read_z(jb,1,jkb) ) * ( zdepth(jk) - pdta_read_z(jb,1,jkb+1) ) <= 0._wp ) &911 & .AND. ( pdta_read_z(jb,1,jkb+1) /= zfv_alt) ) THEN! linear interpolation between 2 levels898 ! 899 ! --- interpolate bdy data on the model grid --- ! 900 DO jk = 1, jpk 901 IF( zdepth(jk) <= pdta_read_z(jb,1,1) ) THEN ! above the first level of external data 902 pdta(jb,1,jk) = pdta_read(jb,1,1) 903 ELSEIF( zdepth(jk) > pdta_read_z(jb,1,ipkmax) ) THEN ! below the last level of external data /= FillValue 904 pdta(jb,1,jk) = pdta_read(jb,1,ipkmax) 905 ELSE ! inbetween: vertical interpolation between jkb & jkb+1 906 DO jkb = 1, ipkmax-1 907 IF( ( ( zdepth(jk) - pdta_read_z(jb,1,jkb) ) * ( zdepth(jk) - pdta_read_z(jb,1,jkb+1) ) ) <= 0._wp ) THEN ! linear interpolation between 2 levels 912 908 zi = ( zdepth(jk) - pdta_read_z(jb,1,jkb) ) / ( pdta_read_z(jb,1,jkb+1) - pdta_read_z(jb,1,jkb) ) 913 pdta(jb,1,jk) = pdta_read(jb,1,jkb) + ( pdta_read (jb,1,jkb+1) - pdta_read (jb,1,jkb) ) * zi909 pdta(jb,1,jk) = pdta_read(jb,1,jkb) + zi * ( pdta_read(jb,1,jkb+1) - pdta_read(jb,1,jkb) ) 914 910 ENDIF 915 911 END DO 916 912 ENDIF 917 END DO ! jpk913 END DO 918 914 ! 919 915 END DO ! ipi 920 921 IF(kgrd == 2) THEN ! do we need to adjust the transport term? 916 917 ! --- mask data and adjust transport --- ! 918 SELECT CASE( kgrd ) 919 920 CASE(1) ! mask data (probably unecessary) 922 921 DO jb = 1, ipi 923 922 ji = idx_bdy(kbdy)%nbi(jb,kgrd) 924 923 jj = idx_bdy(kbdy)%nbj(jb,kgrd) 925 zh = SUM(pdta_read_dz(jb,1,:) ) 924 DO jk = 1, jpk 925 pdta(jb,1,jk) = pdta(jb,1,jk) * tmask(ji,jj,jk) 926 END DO 927 END DO 928 929 CASE(2) ! adjust the U-transport term 930 DO jb = 1, ipi 931 ji = idx_bdy(kbdy)%nbi(jb,kgrd) 932 jj = idx_bdy(kbdy)%nbj(jb,kgrd) 926 933 ztrans = 0._wp 934 DO jkb = 1, ipkb ! calculate transport on input grid 935 IF( pdta_read(jb,1,jkb) /= pfv ) ztrans = ztrans + pdta_read(jb,1,jkb) * pdta_read_dz(jb,1,jkb) 936 ENDDO 927 937 ztrans_new = 0._wp 928 DO jkb = 1, ipkb ! calculate transport on input grid929 ztrans = ztrans + pdta_read(jb,1,jkb) * pdta_read_dz(jb, 1,jkb)930 ENDDO931 938 DO jk = 1, jpk ! calculate transport on model grid 932 ztrans_new = ztrans_new + pdta(jb,1,jk ) * e3u_n(ji,jj,jk) * umask(ji,jj,jk)939 ztrans_new = ztrans_new + pdta(jb,1,jk ) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) 933 940 ENDDO 934 941 DO jk = 1, jpk ! make transport correction 935 942 IF(ldtotvel) THEN ! bdy data are total velocity so adjust bt transport term to match input data 936 943 pdta(jb,1,jk) = ( pdta(jb,1,jk) + ( ztrans - ztrans_new ) * r1_hu_n(ji,jj) ) * umask(ji,jj,jk) 937 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero938 pdta(jb,1,jk) = pdta(jb,1,jk) + ( 0._wp - ztrans_new ) * r1_hu_n(ji,jj)* umask(ji,jj,jk)944 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 945 pdta(jb,1,jk) = ( pdta(jb,1,jk) + ( 0._wp - ztrans_new ) * r1_hu_n(ji,jj) ) * umask(ji,jj,jk) 939 946 ENDIF 940 947 ENDDO 941 948 ENDDO 942 ENDIF 943 944 IF(kgrd == 3) THEN ! do we need to adjust the transport term? 949 950 CASE(3) ! adjust the V-transport term 945 951 DO jb = 1, ipi 946 952 ji = idx_bdy(kbdy)%nbi(jb,kgrd) 947 953 jj = idx_bdy(kbdy)%nbj(jb,kgrd) 948 zh = SUM(pdta_read_dz(jb,1,:) )949 954 ztrans = 0._wp 955 DO jkb = 1, ipkb ! calculate transport on input grid 956 IF( pdta_read(jb,1,jkb) /= pfv ) ztrans = ztrans + pdta_read(jb,1,jkb) * pdta_read_dz(jb,1,jkb) 957 ENDDO 950 958 ztrans_new = 0._wp 951 DO jkb = 1, ipkb ! calculate transport on input grid952 ztrans = ztrans + pdta_read(jb,1,jkb) * pdta_read_dz(jb, 1,jkb)953 ENDDO954 959 DO jk = 1, jpk ! calculate transport on model grid 955 ztrans_new = ztrans_new + pdta(jb,1,jk ) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk)960 ztrans_new = ztrans_new + pdta(jb,1,jk ) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) 956 961 ENDDO 957 962 DO jk = 1, jpk ! make transport correction 958 963 IF(ldtotvel) THEN ! bdy data are total velocity so adjust bt transport term to match input data 959 964 pdta(jb,1,jk) = ( pdta(jb,1,jk) + ( ztrans - ztrans_new ) * r1_hv_n(ji,jj) ) * vmask(ji,jj,jk) 960 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero961 pdta(jb,1,jk) = pdta(jb,1,jk) + ( 0._wp - ztrans_new ) * r1_hv_n(ji,jj)* vmask(ji,jj,jk)965 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 966 pdta(jb,1,jk) = ( pdta(jb,1,jk) + ( 0._wp - ztrans_new ) * r1_hv_n(ji,jj) ) * vmask(ji,jj,jk) 962 967 ENDIF 963 968 ENDDO 964 969 ENDDO 965 END IF966 970 END SELECT 971 967 972 END SUBROUTINE fld_bdy_interp 968 973 969 974 970 975 SUBROUTINE fld_rot( kt, sd ) 971 976 !!---------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
