Ignore:
Timestamp:
2014-06-26T12:22:40+02:00 (7 years ago)
Author:
jamesharle
Message:

Update of fldread to handle depth information in BDY files and addition of an interpolation routine. Updated BDY code to handle T/S BDY interpolation on the fly. Conservative remapping of U/V still to be coded. Not compiled or test yet.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90

    r4371 r4694  
    6767      INTEGER                         ::   nreclast     ! last record to be read in the current file 
    6868      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 
    6971   END TYPE FLD 
    7072 
     
    113115CONTAINS 
    114116 
    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 ) 
    116118      !!--------------------------------------------------------------------- 
    117119      !!                    ***  ROUTINE fld_read  *** 
     
    134136                                                            !   kt_offset = +1 => fields at "after"  time level 
    135137                                                            !   etc. 
     138      !! 
     139      INTEGER  , INTENT(in   ), OPTIONAL     ::   jpk_bdy   ! number of vertical levels in the BDY data 
    136140      !! 
    137141      INTEGER  ::   itmp       ! temporary variable 
     
    169173         DO jf = 1, imf  
    170174            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 
    172180         END DO 
    173181         IF( lwp ) CALL wgt_print()                ! control print 
     
    261269 
    262270               ! 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 
    265276            ENDIF   ! read new data? 
    266277         END DO                                    ! --- end loop over field --- ! 
     
    303314 
    304315 
    305    SUBROUTINE fld_init( kn_fsbc, sdjf, map ) 
     316   SUBROUTINE fld_init( kn_fsbc, sdjf, map , jpk_bdy) 
    306317      !!--------------------------------------------------------------------- 
    307318      !!                    ***  ROUTINE fld_init  *** 
     
    310321      !!               - if time interpolation, read before data  
    311322      !!---------------------------------------------------------------------- 
    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 
    315327      !! 
    316328      LOGICAL :: llprevyr              ! are we reading previous year  file? 
     
    407419 
    408420         ! 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 
    410426 
    411427         clfmt = "('fld_init : time-interpolation for ', a, ' read previous record = ', i6, ' at time = ', f7.2, ' days')" 
     
    577593 
    578594 
    579    SUBROUTINE fld_get( sdjf, map ) 
     595   SUBROUTINE fld_get( sdjf, map, jpk_bdy ) 
    580596      !!--------------------------------------------------------------------- 
    581597      !!                    ***  ROUTINE fld_get  *** 
     
    583599      !! ** Purpose :   read the data 
    584600      !!---------------------------------------------------------------------- 
    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 
    587604      !! 
    588605      INTEGER                  ::   ipk    ! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk ) 
     
    597614      ! 
    598615      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 
    601628         ENDIF 
    602629      ELSE IF( LEN(TRIM(sdjf%wgtname)) > 0 ) THEN 
    603630         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 ) 
    608635         ENDIF 
    609636      ELSE 
     
    654681   END SUBROUTINE fld_get 
    655682 
    656    SUBROUTINE fld_map( num, clvar, dta, nrec, map ) 
     683   SUBROUTINE fld_map( num, clvar, dta, nrec, map, igrd, ibdy, jpk_bdy ) 
    657684      !!--------------------------------------------------------------------- 
    658685      !!                    ***  ROUTINE fld_map  *** 
     
    669696      INTEGER                   , INTENT(in ) ::   nrec    ! record number to read (ie time slice) 
    670697      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 
    671700      !! 
    672701      INTEGER                                 ::   ipi      ! length of boundary data on local process 
     
    677706      INTEGER                                 ::   ib, ik, ji, jj   ! loop counters 
    678707      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 
    680710      !!--------------------------------------------------------------------- 
    681711             
     
    689719#if defined key_bdy 
    690720      ipj = iom_file(num)%dimsz(2,idvar) 
    691       IF (ipj == 1) THEN ! we assume that this is a structured open boundary file 
     721      IF (ipj == 1) THEN  
    692722         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 
    694730         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 
    695737      ENDIF 
    696738#endif 
     
    701743      SELECT CASE( ipk ) 
    702744      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) 
    704774      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 
    707815         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 
    710841            END DO 
    711842         END DO 
     
    714845            jj=1+floor(REAL(map(ib)-1)/REAL(ilendta)) 
    715846            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 
    718874            END DO 
    719875         END DO 
    720876      ENDIF 
    721877 
    722    END SUBROUTINE fld_map 
    723  
     878   END SUBROUTINE fld_bdy_interp 
    724879 
    725880   SUBROUTINE fld_rot( kt, sd ) 
     
    9051060            &   CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdf(jf)%clrootname)//') needs ln_clim = .FALSE.') 
    9061061         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 
    9071064      END DO 
    9081065 
Note: See TracChangeset for help on using the changeset viewer.