Changeset 4694 for branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90
- Timestamp:
- 2014-06-26T12:22:40+02:00 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90
r4371 r4694 67 67 INTEGER :: nreclast ! last record to be read in the current file 68 68 CHARACTER(len = 256) :: lsmname ! current name of the NetCDF mask file acting as a key 69 INTEGER :: igrd ! grid type for bdy data 70 INTEGER :: ibdy ! bdy set id number 69 71 END TYPE FLD 70 72 … … 113 115 CONTAINS 114 116 115 SUBROUTINE fld_read( kt, kn_fsbc, sd, map, kit, kt_offset )117 SUBROUTINE fld_read( kt, kn_fsbc, sd, map, kit, kt_offset, jpk_bdy ) 116 118 !!--------------------------------------------------------------------- 117 119 !! *** ROUTINE fld_read *** … … 134 136 ! kt_offset = +1 => fields at "after" time level 135 137 ! etc. 138 !! 139 INTEGER , INTENT(in ), OPTIONAL :: jpk_bdy ! number of vertical levels in the BDY data 136 140 !! 137 141 INTEGER :: itmp ! temporary variable … … 169 173 DO jf = 1, imf 170 174 IF( PRESENT(map) ) imap = map(jf) 171 CALL fld_init( kn_fsbc, sd(jf), imap ) ! read each before field (put them in after as they will be swapped) 175 IF( PRESENT(jpk_bdy) ) THEN 176 CALL fld_init( kn_fsbc, sd(jf), imap, jpk_bdy ) ! read each before field (put them in after as they will be swapped) 177 ELSE 178 CALL fld_init( kn_fsbc, sd(jf), imap ) ! read each before field (put them in after as they will be swapped) 179 ENDIF 172 180 END DO 173 181 IF( lwp ) CALL wgt_print() ! control print … … 261 269 262 270 ! read after data 263 CALL fld_get( sd(jf), imap ) 264 271 IF( PRESENT(jpk_bdy) ) THEN 272 CALL fld_get( sd(jf), imap, jpk_bdy) 273 ELSE 274 CALL fld_get( sd(jf), imap ) 275 ENDIF 265 276 ENDIF ! read new data? 266 277 END DO ! --- end loop over field --- ! … … 303 314 304 315 305 SUBROUTINE fld_init( kn_fsbc, sdjf, map )316 SUBROUTINE fld_init( kn_fsbc, sdjf, map , jpk_bdy) 306 317 !!--------------------------------------------------------------------- 307 318 !! *** ROUTINE fld_init *** … … 310 321 !! - if time interpolation, read before data 311 322 !!---------------------------------------------------------------------- 312 INTEGER , INTENT(in ) :: kn_fsbc ! sbc computation period (in time step) 313 TYPE(FLD), INTENT(inout) :: sdjf ! input field related variables 314 TYPE(MAP_POINTER),INTENT(in) :: map ! global-to-local mapping indices 323 INTEGER , INTENT(in ) :: kn_fsbc ! sbc computation period (in time step) 324 TYPE(FLD), INTENT(inout) :: sdjf ! input field related variables 325 TYPE(MAP_POINTER),INTENT(in) :: map ! global-to-local mapping indices 326 INTEGER , INTENT(in), OPTIONAL :: jpk_bdy ! number of vertical levels in the BDY data 315 327 !! 316 328 LOGICAL :: llprevyr ! are we reading previous year file? … … 407 419 408 420 ! read before data in after arrays(as we will swap it later) 409 CALL fld_get( sdjf, map ) 421 IF( PRESENT(jpk_bdy) ) THEN 422 CALL fld_get( sdjf, map, jpk_bdy ) 423 ELSE 424 CALL fld_get( sdjf, map ) 425 ENDIF 410 426 411 427 clfmt = "('fld_init : time-interpolation for ', a, ' read previous record = ', i6, ' at time = ', f7.2, ' days')" … … 577 593 578 594 579 SUBROUTINE fld_get( sdjf, map )595 SUBROUTINE fld_get( sdjf, map, jpk_bdy ) 580 596 !!--------------------------------------------------------------------- 581 597 !! *** ROUTINE fld_get *** … … 583 599 !! ** Purpose : read the data 584 600 !!---------------------------------------------------------------------- 585 TYPE(FLD), INTENT(inout) :: sdjf ! input field related variables 586 TYPE(MAP_POINTER),INTENT(in) :: map ! global-to-local mapping indices 601 TYPE(FLD), INTENT(inout) :: sdjf ! input field related variables 602 TYPE(MAP_POINTER),INTENT(in) :: map ! global-to-local mapping indices 603 INTEGER , INTENT(in), OPTIONAL :: jpk_bdy ! number of vertical levels in the bdy data 587 604 !! 588 605 INTEGER :: ipk ! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk ) … … 597 614 ! 598 615 IF( ASSOCIATED(map%ptr) ) THEN 599 IF( sdjf%ln_tint ) THEN ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1), map%ptr ) 600 ELSE ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,: ), sdjf%nrec_a(1), map%ptr ) 616 IF( PRESENT(jpk_bdy) ) THEN 617 IF( sdjf%ln_tint ) THEN ; 618 CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1), map%ptr, sdjf%igrd, sdjf%ibdy, jpk_bdy ) 619 ELSE ; 620 CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,: ), sdjf%nrec_a(1), map%ptr, sdjf%igrd, sdjf%ibdy, jpk_bdy ) 621 ENDIF 622 ELSE 623 IF( sdjf%ln_tint ) THEN ; 624 CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1), map%ptr ) 625 ELSE ; 626 CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,: ), sdjf%nrec_a(1), map%ptr ) 627 ENDIF 601 628 ENDIF 602 629 ELSE IF( LEN(TRIM(sdjf%wgtname)) > 0 ) THEN 603 630 CALL wgt_list( sdjf, iw ) 604 IF( sdjf%ln_tint ) THEN ; CALL fld_interp( sdjf%num, sdjf%clvar, iw , ipk , sdjf%fdta(:,:,:,2), &605 &sdjf%nrec_a(1), sdjf%lsmname )606 ELSE ; CALL fld_interp( sdjf%num, sdjf%clvar, iw , ipk , sdjf%fnow(:,:,: ), &607 &sdjf%nrec_a(1), sdjf%lsmname )631 IF( sdjf%ln_tint ) THEN ; 632 CALL fld_interp( sdjf%num, sdjf%clvar, iw , ipk , sdjf%fdta(:,:,:,2), sdjf%nrec_a(1), sdjf%lsmname ) 633 ELSE ; 634 CALL fld_interp( sdjf%num, sdjf%clvar, iw , ipk , sdjf%fnow(:,:,: ), sdjf%nrec_a(1), sdjf%lsmname ) 608 635 ENDIF 609 636 ELSE … … 654 681 END SUBROUTINE fld_get 655 682 656 SUBROUTINE fld_map( num, clvar, dta, nrec, map )683 SUBROUTINE fld_map( num, clvar, dta, nrec, map, igrd, ibdy, jpk_bdy ) 657 684 !!--------------------------------------------------------------------- 658 685 !! *** ROUTINE fld_map *** … … 669 696 INTEGER , INTENT(in ) :: nrec ! record number to read (ie time slice) 670 697 INTEGER, DIMENSION(:) , INTENT(in ) :: map ! global-to-local mapping indices 698 INTEGER , INTENT(in), OPTIONAL :: igrd, ibdy, jpk_bdy ! grid type, set number and number of vertical levels in the bdy data 699 INTEGER :: jpkm1_bdy! number of vertical levels in the bdy data minus 1 671 700 !! 672 701 INTEGER :: ipi ! length of boundary data on local process … … 677 706 INTEGER :: ib, ik, ji, jj ! loop counters 678 707 INTEGER :: ierr 679 REAL(wp), POINTER, DIMENSION(:,:,:) :: dta_read ! work space for global data 708 REAL(wp), POINTER, DIMENSION(:,:,:) :: dta_read ! work space for global data 709 REAL(wp), POINTER, DIMENSION(:,:,:) :: dta_read_z ! work space for global data 680 710 !!--------------------------------------------------------------------- 681 711 … … 689 719 #if defined key_bdy 690 720 ipj = iom_file(num)%dimsz(2,idvar) 691 IF (ipj == 1) THEN ! we assume that this is a structured open boundary file721 IF (ipj == 1) THEN 692 722 dta_read => dta_global 693 ELSE 723 IF( PRESENT(jpk_bdy) ) THEN 724 IF( jpk_bdy>0 ) THEN 725 dta_read_z => dta_global_z 726 jpkm1_bdy = jpk_bdy-1 727 ENDIF 728 ENDIF 729 ELSE ! we assume that this is a structured open boundary file 694 730 dta_read => dta_global2 731 IF( PRESENT(jpk_bdy) ) THEN 732 IF( jpk_bdy>0 ) THEN 733 dta_read_z => dta_global2_z 734 jpkm1_bdy = jpk_bdy-1 735 ENDIF 736 ENDIF 695 737 ENDIF 696 738 #endif … … 701 743 SELECT CASE( ipk ) 702 744 CASE(1) ; CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1 ), nrec ) 703 CASE DEFAULT ; CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:ipk), nrec ) 745 CASE DEFAULT ; 746 747 IF( PRESENT(jpk_bdy) .AND. jpk_bdy>0 ) THEN ! boundary data not on model grid: vertical interpolation 748 CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:jpk_bdy), nrec ) 749 SELECT CASE( igrd ) 750 CASE(1); CALL iom_get ( num, jpdom_unknown, 'gdept', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) ) 751 CASE(2); CALL iom_get ( num, jpdom_unknown, 'gdepu', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) ) 752 CASE(3); CALL iom_get ( num, jpdom_unknown, 'gdepv', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) ) 753 END SELECT 754 CALL iom_getatt(num, '_FillValue', fv, cdvar=clvar ) 755 CALL fld_bdy_interp(dta_read, dta_read_z, map, jpk_bdy, igrd, ibdy, fv, dta) 756 ELSE ! boundary data assumed to be on model grid 757 CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:ipk), nrec ) 758 IF (ipj==1) THEN ! we assume that this is an un-structured open boundary file 759 DO ib = 1, ipi 760 DO ik = 1, ipk 761 dta(ib,1,ik) = dta_read(map(ib),1,ik) 762 END DO 763 END DO 764 ELSE ! we assume that this is a structured open boundary file 765 DO ib = 1, ipi 766 jj=1+floor(REAL(map(ib)-1)/REAL(ilendta)) 767 ji=map(ib)-(jj-1)*ilendta 768 DO ik = 1, ipk 769 dta(ib,1,ik) = dta_read(ji,jj,ik) 770 END DO 771 END DO 772 ENDIF 773 ENDIF ! PRESENT(jpk_bdy) 704 774 END SELECT 705 ! 706 IF (ipj==1) THEN 775 776 END SUBROUTINE fld_map 777 778 SUBROUTINE fld_bdy_interp(dta_read, dta_read_z, map, jpk_bdy, igrd, ibdy, fv, dta) 779 780 !!--------------------------------------------------------------------- 781 !! *** ROUTINE fld_bdy_interp *** 782 !! 783 !! ** Purpose : on the fly vertical interpolation to allow the use of 784 !! boundary data from non-native vertical grid 785 !!---------------------------------------------------------------------- 786 #if defined key_bdy 787 USE bdy_oce, ONLY: idx_bdy ! indexing for map <-> ij transformation 788 #endif 789 790 REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in ) :: dta_read ! work space for global data 791 REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in ) :: dta_read_z ! work space for global data 792 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: dta ! output field on model grid (2 dimensional) 793 INTEGER, DIMENSION(:) , INTENT(in ) :: map ! global-to-local mapping indices 794 INTEGER , INTENT(in) :: igrd, ib_bdy, jpk_bdy ! number of levels in bdy data 795 INTEGER :: jpkm1_bdy ! number of levels in bdy data minus 1 796 !! 797 INTEGER :: ipi ! length of boundary data on local process 798 INTEGER :: ipj ! length of dummy dimension ( = 1 ) 799 INTEGER :: ipk, ipkm1 ! number of vertical levels of dta ( 2D: ipk=1 ; 3D: ipk=jpk ) 800 INTEGER :: ilendta ! length of data in file 801 INTEGER :: ib, ik, ikk! loop counters 802 REAL(wp) :: zl, zi ! tmp variable for current depth and interpolation factor 803 REAL(wp) :: fv, fv_alt ! fillvalue and alternative -ABS(fv) 804 !!--------------------------------------------------------------------- 805 806 ipi = SIZE( dta, 1 ) 807 ipj = SIZE( dta_read, 2 ) 808 ipk = SIZE( dta, 3 ) 809 ilendta = SIZE( dta_read, 1 ) 810 jpkm1_bdy = jpk_bdy-1 811 812 fv_alt = -ABS(fv) ! set _FillValue < 0 as we make use of MAXVAL and MAXLOC later 813 ! 814 IF (ipj==1) THEN ! we assume that this is an un-structured open boundary file 707 815 DO ib = 1, ipi 708 DO ik = 1, ipk 709 dta(ib,1,ik) = dta_read(map(ib),1,ik) 816 DO ik = 1, ipk 817 IF( ( dta_read(map(ib),1,ik) == fv ) ) THEN 818 dta_read_z(map(ib),1,ik) = fv_alt ! safety: put fillvalue into external depth field so consistent with data 819 ENDIF 820 dta(ib,1,ik) = fv_alt ! put fillvalue into new field as if all goes well all wet points will be replaced 821 ENDDO 822 ENDDO 823 ! 824 DO ib = 1, ipi 825 DO ik = 1, ipk 826 zl = gdept_1(idx_bdy(ib_bdy)%nbi(ib,igrd),idx_bdy(ib_bdy)%nbj(ib,igrd),ik) ! if using in step could use fsdept instead of gdept_1? 827 IF( zl < dta_read_z(map(ib),1,1) ) THEN ! above the first level of external data 828 dta(ib,1,ik) = dta_read(map(ib),1,1) 829 ELSEIF( zl > MAXVAL(dta_read_z(map(ib),1,:),1) ) THEN ! below the last level of external data 830 dta(ib,1,ik) = dta_read(map(ib),1,MAXLOC(dta_read_z(map(ib),1,:),1)) 831 ELSE ! inbetween : vertical interpolation between ikk & ikk+1 832 DO ikk = 1, jpkm1_bdy ! when gdept_1(ikk) < zl < gdept_1(ikk+1) 833 IF( ( (zl-dta_read_z(map(ib),1,ikk)) * (zl-dta_read_z(map(ib),1,ikk+1)) <= 0._wp) & 834 & .AND. (dta_read_z(map(ib),1,ikk+1) /= fv_alt)) THEN 835 zi = ( zl - dta_read_z(map(ib),1,ikk) ) / (dta_read_z(map(ib),1,ikk+1)-dta_read_z(map(ib),1,ikk)) 836 dta(ib,1,ik) = dta_read(map(ib),1,ikk) + & 837 & ( dta_read(map(ib),1,ikk+1) - dta_read(map(ib),1,ikk) ) * zi 838 ENDIF 839 END DO 840 ENDIF 710 841 END DO 711 842 END DO … … 714 845 jj=1+floor(REAL(map(ib)-1)/REAL(ilendta)) 715 846 ji=map(ib)-(jj-1)*ilendta 716 DO ik = 1, ipk 717 dta(ib,1,ik) = dta_read(ji,jj,ik) 847 DO ik = 1, ipk 848 IF( ( dta_read(ji,jj,ik) == fv ) ) THEN 849 dta_read_z(map(ib),1,ik) = fv_alt ! safety: put fillvalue into external depth field so consistent with data 850 ENDIF 851 dta(ib,1,ik) = fv_alt ! put fillvalue into new field as if all goes well all wet points will be replaced 852 ENDDO 853 ENDDO 854 ! 855 DO ib = 1, ipi 856 jj=1+floor(REAL(map(ib)-1)/REAL(ilendta)) 857 ji=map(ib)-(jj-1)*ilendta 858 DO ik = 1, ipk 859 zl = gdept_1(idx_bdy(ib_bdy)%nbi(ib,igrd),idx_bdy(ib_bdy)%nbj(ib,igrd),ik) ! if using in step could use fsdept instead of gdept_1? 860 IF( zl < dta_read_z(ji,jj,1) ) THEN ! above the first level of external data 861 dta(ib,1,ik) = dta_read(ji,jj,1,1) 862 ELSEIF( zl > MAXVAL(dta_read_z(ji,ji,:),1) ) THEN ! below the last level of external data 863 dta(ib,1,ik) = dta_read(ji,jj,MAXLOC(dta_read_z(ji,jj,:),1)) 864 ELSE ! inbetween : vertical interpolation between ikk & ikk+1 865 DO ikk = 1, jpkm1_bdy ! when gdept_1(ikk) < zl < gdept_1(ikk+1) 866 IF( ( (zl-dta_read_z(ji,jj,ikk)) * (zl-dta_read_z(ji,jj,ikk+1)) <= 0._wp) & 867 & .AND. (dta_read_z(ji,jj,ikk+1) /= fv_alt)) THEN 868 zi = ( zl - dta_read_z(ji,jj,ikk) ) / (dta_read_z(ji,jj,ikk+1)-dta_read_z(ji,jj,ikk)) 869 dta(ib,1,ik) = dta_read(ji,jj,ikk) + & 870 & ( dta_read(ji,jj,1,ikk+1) - dta_read(ji,jj,ikk) ) * zi 871 ENDIF 872 END DO 873 ENDIF 718 874 END DO 719 875 END DO 720 876 ENDIF 721 877 722 END SUBROUTINE fld_map 723 878 END SUBROUTINE fld_bdy_interp 724 879 725 880 SUBROUTINE fld_rot( kt, sd ) … … 905 1060 & CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdf(jf)%clrootname)//') needs ln_clim = .FALSE.') 906 1061 sdf(jf)%nreclast = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during fld_clopn 1062 sdf(jf)%igrd = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during bdy_dta_init 1063 sdf(jf)%ibdy = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during bdy_dta_init 907 1064 END DO 908 1065
Note: See TracChangeset
for help on using the changeset viewer.
