Changeset 12367


Ignore:
Timestamp:
2020-02-11T19:37:34+01:00 (8 weeks ago)
Author:
clem
Message:

solve ticket #2380

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/src/OCE/SBC/fldread.F90

    r12248 r12367  
    833833 
    834834      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdta_read       ! data read in bdy file 
    835       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdta_read_z     ! depth of the data read in bdy file 
    836       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdta_read_dz    ! thickness of the levels in bdy file 
     835      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 
    837837      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pdta            ! output field on model grid (2 dimensional) 
    838838      REAL(wp)                  , INTENT(in   ) ::   pfv             ! fillvalue of the data read in bdy file 
     
    841841      INTEGER                   , INTENT(in   ) ::   kbdy            ! bdy number 
    842842      !! 
    843       INTEGER                                   ::   ipi             ! length of boundary data on local process 
    844       INTEGER                                   ::   ipkb            ! number of vertical levels in boundary data file 
    845       INTEGER                                   ::   jb, ji, jj, jk, jkb   ! loop counters 
    846       REAL(wp)                                  ::   zcoef 
    847       REAL(wp)                                  ::   zl, zi, zh      ! tmp variable for current depth and interpolation factor 
    848       REAL(wp)                                  ::   zfv_alt, ztrans, ztrans_new ! fillvalue and alternative -ABS(pfv) 
    849       REAL(wp), DIMENSION(jpk)                  ::   zdepth, zdhalf  ! level and half-level depth 
     843      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 
    850850      !!--------------------------------------------------------------------- 
    851851      
     
    853853      ipkb = SIZE( pdta_read, 3 ) 
    854854       
    855       zfv_alt = -ABS(pfv)  ! set _FillValue < 0 as we make use of MAXVAL and MAXLOC later 
    856       ! 
    857       WHERE( pdta_read == pfv ) 
    858          pdta_read_z  = zfv_alt ! safety: put fillvalue into external depth field so consistent with data 
    859          pdta_read_dz = 0._wp   ! safety: put 0._wp into external thickness factors to ensure transport is correct 
    860       ENDWHERE 
    861        
    862855      DO jb = 1, ipi 
    863856         ji = idx_bdy(kbdy)%nbi(jb,kgrd) 
    864857         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 --- ! 
    867866         SELECT CASE( kgrd )                          
    868          CASE(1) 
    869             ! depth of T points: 
     867         CASE(1)            ! depth of T points: 
    870868            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... 
    874871            zdhalf(1) = 0.0_wp 
    875872            zdepth(1) = 0.5_wp * e3uw_n(ji,jj,1) 
     
    881878               zcoef = ( umask(ji,jj,jk) - wumask(ji,jj,jk) ) 
    882879               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)) 
    885882            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... 
    889885            zdhalf(1) = 0.0_wp 
    890886            zdepth(1) = 0.5_wp * e3vw_n(ji,jj,1) 
     
    896892               zcoef = ( vmask(ji,jj,jk) - wvmask(ji,jj,jk) ) 
    897893               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)) 
    900896            END DO 
    901897         END SELECT 
    902          ! 
    903          DO jk = 1, jpk                       
    904             IF(     zdepth(jk) < pdta_read_z(jb,1,          1) ) THEN                ! above the first level of external data 
    905                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 data  
    907                pdta(jb,1,jk) =  pdta_read(jb,1,MAXLOC(pdta_read_z(jb,1,:),1)) 
    908             ELSE                                                             ! inbetween: vertical interpolation between jkb & jkb+1 
    909                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 levels 
     898         !          
     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 
    912908                     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) ) * zi 
     909                     pdta(jb,1,jk) = pdta_read(jb,1,jkb) + zi * ( pdta_read(jb,1,jkb+1) - pdta_read(jb,1,jkb) ) 
    914910                  ENDIF 
    915911               END DO 
    916912            ENDIF 
    917          END DO   ! jpk 
     913         END DO 
    918914         ! 
    919915      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) 
    922921         DO jb = 1, ipi 
    923922            ji = idx_bdy(kbdy)%nbi(jb,kgrd) 
    924923            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) 
    926933            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 
    927937            ztrans_new = 0._wp 
    928             DO jkb = 1, ipkb                              ! calculate transport on input grid 
    929                ztrans     = ztrans     + pdta_read(jb,1,jkb) * pdta_read_dz(jb, 1,jkb) 
    930             ENDDO 
    931938            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) 
    933940            ENDDO 
    934941            DO jk = 1, jpk                                ! make transport correction 
    935942               IF(ldtotvel) THEN ! bdy data are total velocity so adjust bt transport term to match input data 
    936943                  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 zero 
    938                   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) 
    939946               ENDIF 
    940947            ENDDO 
    941948         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 
    945951         DO jb = 1, ipi 
    946952            ji = idx_bdy(kbdy)%nbi(jb,kgrd) 
    947953            jj = idx_bdy(kbdy)%nbj(jb,kgrd) 
    948             zh  = SUM(pdta_read_dz(jb,1,:) ) 
    949954            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 
    950958            ztrans_new = 0._wp 
    951             DO jkb = 1, ipkb                              ! calculate transport on input grid 
    952                ztrans     = ztrans     + pdta_read(jb,1,jkb) * pdta_read_dz(jb, 1,jkb) 
    953             ENDDO 
    954959            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) 
    956961            ENDDO 
    957962            DO jk = 1, jpk                                ! make transport correction 
    958963               IF(ldtotvel) THEN ! bdy data are total velocity so adjust bt transport term to match input data 
    959964                  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 zero 
    961                   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) 
    962967               ENDIF 
    963968            ENDDO 
    964969         ENDDO 
    965       ENDIF 
    966        
     970      END SELECT 
     971 
    967972   END SUBROUTINE fld_bdy_interp 
    968973 
    969  
     974    
    970975   SUBROUTINE fld_rot( kt, sd ) 
    971976      !!--------------------------------------------------------------------- 
Note: See TracChangeset for help on using the changeset viewer.