Changeset 11480 for NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/fldread.F90
- Timestamp:
- 2019-08-29T11:23:25+02:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/fldread.F90
r10922 r11480 130 130 CONTAINS 131 131 132 SUBROUTINE fld_read( kt, kn_fsbc, sd, map, kit, kt_offset, jpk_bdy, fvl )132 SUBROUTINE fld_read( kt, kn_fsbc, sd, map, kit, kt_offset, jpk_bdy, fvl, Kmm ) 133 133 !!--------------------------------------------------------------------- 134 134 !! *** ROUTINE fld_read *** … … 153 153 INTEGER , INTENT(in ), OPTIONAL :: jpk_bdy ! number of vertical levels in the BDY data 154 154 LOGICAL , INTENT(in ), OPTIONAL :: fvl ! number of vertical levels in the BDY data 155 INTEGER , INTENT(in ), OPTIONAL :: Kmm ! ocean time level index 155 156 !! 156 157 INTEGER :: itmp ! local variable … … 287 288 ! read after data 288 289 IF( PRESENT(jpk_bdy) ) THEN 289 CALL fld_get( sd(jf), imap, jpk_bdy, fvl )290 CALL fld_get( sd(jf), imap, jpk_bdy, fvl, Kmm ) 290 291 ELSE 291 292 CALL fld_get( sd(jf), imap ) … … 614 615 615 616 616 SUBROUTINE fld_get( sdjf, map, jpk_bdy, fvl )617 SUBROUTINE fld_get( sdjf, map, jpk_bdy, fvl, Kmm ) 617 618 !!--------------------------------------------------------------------- 618 619 !! *** ROUTINE fld_get *** … … 624 625 INTEGER , INTENT(in), OPTIONAL :: jpk_bdy ! number of vertical levels in the bdy data 625 626 LOGICAL , INTENT(in), OPTIONAL :: fvl ! number of vertical levels in the bdy data 627 INTEGER , INTENT(in), OPTIONAL :: Kmm ! ocean time level index 626 628 ! 627 629 INTEGER :: ipk ! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk ) … … 638 640 IF( PRESENT(jpk_bdy) ) THEN 639 641 IF( sdjf%ln_tint ) THEN ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2), & 640 sdjf%nrec_a(1), map, sdjf%igrd, sdjf%ibdy, jpk_bdy, fvl )642 sdjf%nrec_a(1), map, sdjf%igrd, sdjf%ibdy, jpk_bdy, fvl, Kmm ) 641 643 ELSE ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,: ), & 642 sdjf%nrec_a(1), map, sdjf%igrd, sdjf%ibdy, jpk_bdy, fvl )644 sdjf%nrec_a(1), map, sdjf%igrd, sdjf%ibdy, jpk_bdy, fvl, Kmm ) 643 645 ENDIF 644 646 ELSE … … 701 703 END SUBROUTINE fld_get 702 704 703 SUBROUTINE fld_map( num, clvar, dta, nrec, map, igrd, ibdy, jpk_bdy, fvl )705 SUBROUTINE fld_map( num, clvar, dta, nrec, map, igrd, ibdy, jpk_bdy, fvl, Kmm ) 704 706 !!--------------------------------------------------------------------- 705 707 !! *** ROUTINE fld_map *** … … 718 720 INTEGER , INTENT(in), OPTIONAL :: igrd, ibdy, jpk_bdy ! grid type, set number and number of vertical levels in the bdy data 719 721 LOGICAL , INTENT(in), OPTIONAL :: fvl ! grid type, set number and number of vertical levels in the bdy data 722 INTEGER , INTENT(in), OPTIONAL :: Kmm ! ocean time level index 720 723 INTEGER :: jpkm1_bdy! number of vertical levels in the bdy data minus 1 721 724 !! … … 813 816 814 817 IF ( ln_bdy ) & 815 CALL fld_bdy_interp(dta_read, dta_read_z, dta_read_dz, map, jpk_bdy, igrd, ibdy, fv, dta, fvl, ilendta )818 CALL fld_bdy_interp(dta_read, dta_read_z, dta_read_dz, map, jpk_bdy, igrd, ibdy, fv, dta, fvl, ilendta, Kmm) 816 819 817 820 ELSE ! boundary data assumed to be on model grid … … 838 841 END SUBROUTINE fld_map 839 842 840 SUBROUTINE fld_bdy_interp(dta_read, dta_read_z, dta_read_dz, map, jpk_bdy, igrd, ibdy, fv, dta, fvl, ilendta )843 SUBROUTINE fld_bdy_interp(dta_read, dta_read_z, dta_read_dz, map, jpk_bdy, igrd, ibdy, fv, dta, fvl, ilendta, Kmm) 841 844 842 845 !!--------------------------------------------------------------------- … … 857 860 INTEGER , INTENT(in) :: igrd, ibdy, jpk_bdy ! number of levels in bdy data 858 861 INTEGER , INTENT(in) :: ilendta ! length of data in file 862 INTEGER , INTENT(in), OPTIONAL :: Kmm ! ocean time level index 859 863 !! 860 864 INTEGER :: ipi ! length of boundary data on local process … … 900 904 SELECT CASE( igrd ) 901 905 CASE(1) 902 IF( ABS( (zh - ht _n(zij,zjj)) / ht_n(zij,zjj)) * tmask(zij,zjj,1) > 0.01_wp ) THEN906 IF( ABS( (zh - ht(zij,zjj)) / ht(zij,zjj)) * tmask(zij,zjj,1) > 0.01_wp ) THEN 903 907 WRITE(ibstr,"(I10.10)") map%ptr(ib) 904 908 CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 905 ! IF(lwp) WRITE(*,*) 'DEPTHT', zh, sum(e3t(zij,zjj,:, nfld_Nnn), mask=tmask(zij,zjj,:)==1), ht_n(zij,zjj), map%ptr(ib), ib, zij, zjj909 ! IF(lwp) WRITE(*,*) 'DEPTHT', zh, sum(e3t(zij,zjj,:,Kmm), mask=tmask(zij,zjj,:)==1), ht(zij,zjj), map%ptr(ib), ib, zij, zjj 906 910 ENDIF 907 911 CASE(2) 908 IF( ABS( (zh - hu _n(zij,zjj)) * r1_hu_n(zij,zjj)) * umask(zij,zjj,1) > 0.01_wp ) THEN912 IF( ABS( (zh - hu(zij,zjj,Kmm)) * r1_hu(zij,zjj,Kmm)) * umask(zij,zjj,1) > 0.01_wp ) THEN 909 913 WRITE(ibstr,"(I10.10)") map%ptr(ib) 910 914 CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 911 IF(lwp) WRITE(*,*) 'DEPTHU', zh, sum(e3u(zij,zjj,:, nfld_Nnn), mask=umask(zij,zjj,:)==1), sum(umask(zij,zjj,:)), &912 & hu _n(zij,zjj), map%ptr(ib), ib, zij, zjj, narea-1 , &915 IF(lwp) WRITE(*,*) 'DEPTHU', zh, sum(e3u(zij,zjj,:,Kmm), mask=umask(zij,zjj,:)==1), sum(umask(zij,zjj,:)), & 916 & hu(zij,zjj,Kmm), map%ptr(ib), ib, zij, zjj, narea-1 , & 913 917 & dta_read(map%ptr(ib),1,:) 914 918 ENDIF 915 919 CASE(3) 916 IF( ABS( (zh - hv _n(zij,zjj)) * r1_hv_n(zij,zjj)) * vmask(zij,zjj,1) > 0.01_wp ) THEN920 IF( ABS( (zh - hv(zij,zjj,Kmm)) * r1_hv(zij,zjj,Kmm)) * vmask(zij,zjj,1) > 0.01_wp ) THEN 917 921 WRITE(ibstr,"(I10.10)") map%ptr(ib) 918 922 CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') … … 922 926 SELECT CASE( igrd ) 923 927 CASE(1) 924 zl = gdept(zij,zjj,ik, nfld_Nnn) ! if using in step could use fsdept instead of gdept_n?928 zl = gdept(zij,zjj,ik,Kmm) ! if using in step could use fsdept instead of gdept_n? 925 929 CASE(2) 926 930 IF(ln_sco) THEN 927 zl = ( gdept(zij,zjj,ik, nfld_Nnn) + gdept(zij+1,zjj,ik,nfld_Nnn) ) * 0.5_wp ! if using in step could use fsdept instead of gdept_n?931 zl = ( gdept(zij,zjj,ik,Kmm) + gdept(zij+1,zjj,ik,Kmm) ) * 0.5_wp ! if using in step could use fsdept instead of gdept_n? 928 932 ELSE 929 zl = MIN( gdept(zij,zjj,ik, nfld_Nnn), gdept(zij+1,zjj,ik,nfld_Nnn) )933 zl = MIN( gdept(zij,zjj,ik,Kmm), gdept(zij+1,zjj,ik,Kmm) ) 930 934 ENDIF 931 935 CASE(3) 932 936 IF(ln_sco) THEN 933 zl = ( gdept(zij,zjj,ik, nfld_Nnn) + gdept(zij,zjj+1,ik,nfld_Nnn) ) * 0.5_wp ! if using in step could use fsdept instead of gdept_n?937 zl = ( gdept(zij,zjj,ik,Kmm) + gdept(zij,zjj+1,ik,Kmm) ) * 0.5_wp ! if using in step could use fsdept instead of gdept_n? 934 938 ELSE 935 zl = MIN( gdept(zij,zjj,ik, nfld_Nnn), gdept(zij,zjj+1,ik,nfld_Nnn) )939 zl = MIN( gdept(zij,zjj,ik,Kmm), gdept(zij,zjj+1,ik,Kmm) ) 936 940 ENDIF 937 941 END SELECT … … 941 945 dta(ib,1,ik) = dta_read(map%ptr(ib),1,MAXLOC(dta_read_z(map%ptr(ib),1,:),1)) 942 946 ELSE ! inbetween : vertical interpolation between ikk & ikk+1 943 DO ikk = 1, jpkm1_bdy ! when gdept(ikk, nfld_Nnn) < zl < gdept(ikk+1,nfld_Nnn)947 DO ikk = 1, jpkm1_bdy ! when gdept(ikk,Kmm) < zl < gdept(ikk+1,Kmm) 944 948 IF( ( (zl-dta_read_z(map%ptr(ib),1,ikk)) * (zl-dta_read_z(map%ptr(ib),1,ikk+1)) <= 0._wp) & 945 949 & .AND. (dta_read_z(map%ptr(ib),1,ikk+1) /= fv_alt)) THEN … … 965 969 ENDDO 966 970 DO ik = 1, ipk ! calculate transport on model grid 967 ztrans_new = ztrans_new + dta(ib,1,ik) * e3u(zij,zjj,ik, nfld_Nnn) * umask(zij,zjj,ik)971 ztrans_new = ztrans_new + dta(ib,1,ik) * e3u(zij,zjj,ik,Kmm) * umask(zij,zjj,ik) 968 972 ENDDO 969 973 DO ik = 1, ipk ! make transport correction 970 974 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 971 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hu _n(zij,zjj) ) * umask(zij,zjj,ik)975 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hu(zij,zjj,Kmm) ) * umask(zij,zjj,ik) 972 976 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 973 IF( ABS(ztrans * r1_hu _n(zij,zjj)) > 0.01_wp ) &977 IF( ABS(ztrans * r1_hu(zij,zjj,Kmm)) > 0.01_wp ) & 974 978 & CALL ctl_warn('fld_bdy_interp: barotropic component of > 0.01 ms-1 found in baroclinic velocities at') 975 dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hu _n(zij,zjj) * umask(zij,zjj,ik)979 dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hu(zij,zjj,Kmm) * umask(zij,zjj,ik) 976 980 ENDIF 977 981 ENDDO … … 990 994 ENDDO 991 995 DO ik = 1, ipk ! calculate transport on model grid 992 ztrans_new = ztrans_new + dta(ib,1,ik) * e3v(zij,zjj,ik, nfld_Nnn) * vmask(zij,zjj,ik)996 ztrans_new = ztrans_new + dta(ib,1,ik) * e3v(zij,zjj,ik,Kmm) * vmask(zij,zjj,ik) 993 997 ENDDO 994 998 DO ik = 1, ipk ! make transport correction 995 999 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 996 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hv _n(zij,zjj) ) * vmask(zij,zjj,ik)1000 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hv(zij,zjj,Kmm) ) * vmask(zij,zjj,ik) 997 1001 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 998 dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hv _n(zij,zjj) * vmask(zij,zjj,ik)1002 dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hv(zij,zjj,Kmm) * vmask(zij,zjj,ik) 999 1003 ENDIF 1000 1004 ENDDO … … 1025 1029 SELECT CASE( igrd ) 1026 1030 CASE(1) 1027 IF( ABS( (zh - ht _n(zij,zjj)) / ht_n(zij,zjj)) * tmask(zij,zjj,1) > 0.01_wp ) THEN1031 IF( ABS( (zh - ht(zij,zjj)) / ht(zij,zjj)) * tmask(zij,zjj,1) > 0.01_wp ) THEN 1028 1032 WRITE(ibstr,"(I10.10)") map%ptr(ib) 1029 1033 CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 1030 ! IF(lwp) WRITE(*,*) 'DEPTHT', zh, sum(e3t(zij,zjj,:, nfld_Nnn), mask=tmask(zij,zjj,:)==1), ht_n(zij,zjj), map%ptr(ib), ib, zij, zjj1034 ! IF(lwp) WRITE(*,*) 'DEPTHT', zh, sum(e3t(zij,zjj,:,Kmm), mask=tmask(zij,zjj,:)==1), ht(zij,zjj), map%ptr(ib), ib, zij, zjj 1031 1035 ENDIF 1032 1036 CASE(2) 1033 IF( ABS( (zh - hu _n(zij,zjj)) * r1_hu_n(zij,zjj)) * umask(zij,zjj,1) > 0.01_wp ) THEN1037 IF( ABS( (zh - hu(zij,zjj,Kmm)) * r1_hu(zij,zjj,Kmm)) * umask(zij,zjj,1) > 0.01_wp ) THEN 1034 1038 WRITE(ibstr,"(I10.10)") map%ptr(ib) 1035 1039 CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 1036 1040 ENDIF 1037 1041 CASE(3) 1038 IF( ABS( (zh - hv _n(zij,zjj)) * r1_hv_n(zij,zjj)) * vmask(zij,zjj,1) > 0.01_wp ) THEN1042 IF( ABS( (zh - hv(zij,zjj,Kmm)) * r1_hv(zij,zjj,Kmm)) * vmask(zij,zjj,1) > 0.01_wp ) THEN 1039 1043 WRITE(ibstr,"(I10.10)") map%ptr(ib) 1040 1044 CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') … … 1044 1048 SELECT CASE( igrd ) ! coded for sco - need zco and zps option using min 1045 1049 CASE(1) 1046 zl = gdept(zij,zjj,ik, nfld_Nnn) ! if using in step could use fsdept instead of gdept_n?1050 zl = gdept(zij,zjj,ik,Kmm) ! if using in step could use fsdept instead of gdept_n? 1047 1051 CASE(2) 1048 1052 IF(ln_sco) THEN 1049 zl = ( gdept(zij,zjj,ik, nfld_Nnn) + gdept(zij+1,zjj,ik,nfld_Nnn) ) * 0.5_wp ! if using in step could use fsdept instead of gdept_n?1053 zl = ( gdept(zij,zjj,ik,Kmm) + gdept(zij+1,zjj,ik,Kmm) ) * 0.5_wp ! if using in step could use fsdept instead of gdept_n? 1050 1054 ELSE 1051 zl = MIN( gdept(zij,zjj,ik, nfld_Nnn), gdept(zij+1,zjj,ik,nfld_Nnn) )1055 zl = MIN( gdept(zij,zjj,ik,Kmm), gdept(zij+1,zjj,ik,Kmm) ) 1052 1056 ENDIF 1053 1057 CASE(3) 1054 1058 IF(ln_sco) THEN 1055 zl = ( gdept(zij,zjj,ik, nfld_Nnn) + gdept(zij,zjj+1,ik,nfld_Nnn) ) * 0.5_wp ! if using in step could use fsdept instead of gdept_n?1059 zl = ( gdept(zij,zjj,ik,Kmm) + gdept(zij,zjj+1,ik,Kmm) ) * 0.5_wp ! if using in step could use fsdept instead of gdept_n? 1056 1060 ELSE 1057 zl = MIN( gdept(zij,zjj,ik, nfld_Nnn), gdept(zij,zjj+1,ik,nfld_Nnn) )1061 zl = MIN( gdept(zij,zjj,ik,Kmm), gdept(zij,zjj+1,ik,Kmm) ) 1058 1062 ENDIF 1059 1063 END SELECT … … 1063 1067 dta(ib,1,ik) = dta_read(ji,jj,MAXLOC(dta_read_z(ji,jj,:),1)) 1064 1068 ELSE ! inbetween : vertical interpolation between ikk & ikk+1 1065 DO ikk = 1, jpkm1_bdy ! when gdept(ikk, nfld_Nnn) < zl < gdept(ikk+1,nfld_Nnn)1069 DO ikk = 1, jpkm1_bdy ! when gdept(ikk,Kmm) < zl < gdept(ikk+1,Kmm) 1066 1070 IF( ( (zl-dta_read_z(ji,jj,ikk)) * (zl-dta_read_z(ji,jj,ikk+1)) <= 0._wp) & 1067 1071 & .AND. (dta_read_z(ji,jj,ikk+1) /= fv_alt)) THEN … … 1089 1093 ENDDO 1090 1094 DO ik = 1, ipk ! calculate transport on model grid 1091 ztrans_new = ztrans_new + dta(ib,1,ik) * e3u(zij,zjj,ik, nfld_Nnn) * umask(zij,zjj,ik)1095 ztrans_new = ztrans_new + dta(ib,1,ik) * e3u(zij,zjj,ik,Kmm) * umask(zij,zjj,ik) 1092 1096 ENDDO 1093 1097 DO ik = 1, ipk ! make transport correction 1094 1098 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 1095 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hu _n(zij,zjj) ) * umask(zij,zjj,ik)1099 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hu(zij,zjj,Kmm) ) * umask(zij,zjj,ik) 1096 1100 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 1097 dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hu _n(zij,zjj) ) * umask(zij,zjj,ik)1101 dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hu(zij,zjj,Kmm) ) * umask(zij,zjj,ik) 1098 1102 ENDIF 1099 1103 ENDDO … … 1114 1118 ENDDO 1115 1119 DO ik = 1, ipk ! calculate transport on model grid 1116 ztrans_new = ztrans_new + dta(ib,1,ik) * e3v(zij,zjj,ik, nfld_Nnn) * vmask(zij,zjj,ik)1120 ztrans_new = ztrans_new + dta(ib,1,ik) * e3v(zij,zjj,ik,Kmm) * vmask(zij,zjj,ik) 1117 1121 ENDDO 1118 1122 DO ik = 1, ipk ! make transport correction 1119 1123 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 1120 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hv _n(zij,zjj) ) * vmask(zij,zjj,ik)1124 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hv(zij,zjj,Kmm) ) * vmask(zij,zjj,ik) 1121 1125 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 1122 dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hv _n(zij,zjj) ) * vmask(zij,zjj,ik)1126 dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hv(zij,zjj,Kmm) ) * vmask(zij,zjj,ik) 1123 1127 ENDIF 1124 1128 ENDDO
Note: See TracChangeset
for help on using the changeset viewer.