- Timestamp:
- 2020-02-12T13:37:21+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/SBC/fldread.F90
r12350 r12372 571 571 572 572 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pdta_read ! data read in bdy file 573 REAL(wp), DIMENSION(:,:,:), INTENT(in out) :: pdta_read_z ! depth of the data read in bdy file574 REAL(wp), DIMENSION(:,:,:), INTENT(in out) :: pdta_read_dz ! thickness of the levels in bdy file573 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pdta_read_z ! depth of the data read in bdy file 574 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pdta_read_dz ! thickness of the levels in bdy file 575 575 REAL(wp), DIMENSION(:,:,:), INTENT( out) :: pdta ! output field on model grid (2 dimensional) 576 576 REAL(wp) , INTENT(in ) :: pfv ! fillvalue of the data read in bdy file … … 580 580 INTEGER, OPTIONAL , INTENT(in ) :: Kmm ! ocean time level index 581 581 !! 582 INTEGER :: ipi! length of boundary data on local process583 INTEGER :: ipkb! number of vertical levels in boundary data file584 INTEGER :: jb, ji, jj, jk, jkb ! loop counters585 REAL(wp) :: zcoef586 REAL(wp) :: zl, zi, zh ! tmp variable for current depth and interpolation factor587 REAL(wp) :: zfv_alt, ztrans, ztrans_new ! fillvalue and alternative -ABS(pfv)588 REAL(wp), DIMENSION(jpk) :: zdepth, zdhalf! level and half-level depth582 INTEGER :: ipi ! length of boundary data on local process 583 INTEGER :: ipkb ! number of vertical levels in boundary data file 584 INTEGER :: ipkmax ! number of vertical levels in boundary data file where no mask 585 INTEGER :: jb, ji, jj, jk, jkb ! loop counters 586 REAL(wp) :: zcoef, zi ! 587 REAL(wp) :: ztrans, ztrans_new ! transports 588 REAL(wp), DIMENSION(jpk) :: zdepth, zdhalf ! level and half-level depth 589 589 !!--------------------------------------------------------------------- 590 590 … … 592 592 ipkb = SIZE( pdta_read, 3 ) 593 593 594 zfv_alt = -ABS(pfv) ! set _FillValue < 0 as we make use of MAXVAL and MAXLOC later595 !596 WHERE( pdta_read == pfv )597 pdta_read_z = zfv_alt ! safety: put fillvalue into external depth field so consistent with data598 pdta_read_dz = 0._wp ! safety: put 0._wp into external thickness factors to ensure transport is correct599 ENDWHERE600 601 594 DO jb = 1, ipi 602 595 ji = idx_bdy(kbdy)%nbi(jb,kgrd) 603 596 jj = idx_bdy(kbdy)%nbj(jb,kgrd) 604 zh = SUM(pdta_read_dz(jb,1,:) ) 605 ! 597 ! 598 ! --- max jk where input data /= FillValue --- ! 599 ipkmax = 1 600 DO jkb = 2, ipkb 601 IF( pdta_read(jb,1,jkb) /= pfv ) ipkmax = MAX( ipkmax, jkb ) 602 END DO 603 ! 604 ! --- calculate depth at t,u,v points --- ! 606 605 SELECT CASE( kgrd ) 607 CASE(1) 608 ! depth of T points: 606 CASE(1) ! depth of T points: 609 607 zdepth(:) = gdept(ji,jj,:,Kmm) 610 CASE(2) 611 ! depth of U points: we must not use gdept_n as we don't want to do a communication 612 ! --> copy what is done for gdept_n in domvvl... 608 CASE(2) ! depth of U points: we must not use gdept_n as we don't want to do a communication 609 ! --> copy what is done for gdept_n in domvvl... 613 610 zdhalf(1) = 0.0_wp 614 611 zdepth(1) = 0.5_wp * e3uw(ji,jj,1,Kmm) … … 620 617 zcoef = ( umask(ji,jj,jk) - wumask(ji,jj,jk) ) 621 618 zdhalf(jk) = zdhalf(jk-1) + e3u(ji,jj,jk-1,Kmm) 622 zdepth(jk) = zcoef * ( zdhalf(jk ) + 0.5* e3uw(ji,jj,jk,Kmm)) &623 & + (1 -zcoef) * ( zdepth(jk-1) +e3uw(ji,jj,jk,Kmm))619 zdepth(jk) = zcoef * ( zdhalf(jk ) + 0.5_wp * e3uw(ji,jj,jk,Kmm)) & 620 & + (1._wp-zcoef) * ( zdepth(jk-1) + e3uw(ji,jj,jk,Kmm)) 624 621 END DO 625 CASE(3) 626 ! depth of V points: we must not use gdept_n as we don't want to do a communication 627 ! --> copy what is done for gdept_n in domvvl... 622 CASE(3) ! depth of V points: we must not use gdept_n as we don't want to do a communication 623 ! --> copy what is done for gdept_n in domvvl... 628 624 zdhalf(1) = 0.0_wp 629 625 zdepth(1) = 0.5_wp * e3vw(ji,jj,1,Kmm) … … 635 631 zcoef = ( vmask(ji,jj,jk) - wvmask(ji,jj,jk) ) 636 632 zdhalf(jk) = zdhalf(jk-1) + e3v(ji,jj,jk-1,Kmm) 637 zdepth(jk) = zcoef * ( zdhalf(jk ) + 0.5* e3vw(ji,jj,jk,Kmm)) &638 & + (1 -zcoef) * ( zdepth(jk-1) +e3vw(ji,jj,jk,Kmm))633 zdepth(jk) = zcoef * ( zdhalf(jk ) + 0.5_wp * e3vw(ji,jj,jk,Kmm)) & 634 & + (1._wp-zcoef) * ( zdepth(jk-1) + e3vw(ji,jj,jk,Kmm)) 639 635 END DO 640 636 END SELECT 641 ! 642 DO jk = 1, jpk643 IF( zdepth(jk) < pdta_read_z(jb,1, 1) ) THEN ! above the first level of external data644 pdta(jb,1,jk) = pdta_read(jb,1,1)645 ELSEIF( zdepth(jk) > pdta_read_z(jb,1,ipkb) ) THEN ! below the last level of external data646 pdta(jb,1,jk) = pdta_read(jb,1,MAXLOC(pdta_read_z(jb,1,:),1))647 ELSE ! inbetween: vertical interpolation between jkb & jkb+1648 DO jkb = 1, ipkb-1 ! when gdept_n(jkb) < zdepth(jk) < gdept_n(jkb+1)649 IF( ( ( zdepth(jk) - pdta_read_z(jb,1,jkb) ) * ( zdepth(jk) - pdta_read_z(jb,1,jkb+1) ) <= 0._wp ) &650 & .AND. ( pdta_read_z(jb,1,jkb+1) /= zfv_alt) ) THEN! linear interpolation between 2 levels637 ! 638 ! --- interpolate bdy data on the model grid --- ! 639 DO jk = 1, jpk 640 IF( zdepth(jk) <= pdta_read_z(jb,1,1) ) THEN ! above the first level of external data 641 pdta(jb,1,jk) = pdta_read(jb,1,1) 642 ELSEIF( zdepth(jk) > pdta_read_z(jb,1,ipkmax) ) THEN ! below the last level of external data /= FillValue 643 pdta(jb,1,jk) = pdta_read(jb,1,ipkmax) 644 ELSE ! inbetween: vertical interpolation between jkb & jkb+1 645 DO jkb = 1, ipkmax-1 646 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 651 647 zi = ( zdepth(jk) - pdta_read_z(jb,1,jkb) ) / ( pdta_read_z(jb,1,jkb+1) - pdta_read_z(jb,1,jkb) ) 652 pdta(jb,1,jk) = pdta_read(jb,1,jkb) + ( pdta_read (jb,1,jkb+1) - pdta_read (jb,1,jkb) ) * zi648 pdta(jb,1,jk) = pdta_read(jb,1,jkb) + zi * ( pdta_read(jb,1,jkb+1) - pdta_read(jb,1,jkb) ) 653 649 ENDIF 654 650 END DO … … 657 653 ! 658 654 END DO ! ipi 659 660 IF(kgrd == 2) THEN ! do we need to adjust the transport term? 655 656 ! --- mask data and adjust transport --- ! 657 SELECT CASE( kgrd ) 658 659 CASE(1) ! mask data (probably unecessary) 661 660 DO jb = 1, ipi 662 661 ji = idx_bdy(kbdy)%nbi(jb,kgrd) 663 662 jj = idx_bdy(kbdy)%nbj(jb,kgrd) 664 zh = SUM(pdta_read_dz(jb,1,:) ) 663 DO jk = 1, jpk 664 pdta(jb,1,jk) = pdta(jb,1,jk) * tmask(ji,jj,jk) 665 END DO 666 END DO 667 668 CASE(2) ! adjust the U-transport term 669 DO jb = 1, ipi 670 ji = idx_bdy(kbdy)%nbi(jb,kgrd) 671 jj = idx_bdy(kbdy)%nbj(jb,kgrd) 665 672 ztrans = 0._wp 673 DO jkb = 1, ipkb ! calculate transport on input grid 674 IF( pdta_read(jb,1,jkb) /= pfv ) ztrans = ztrans + pdta_read(jb,1,jkb) * pdta_read_dz(jb,1,jkb) 675 ENDDO 666 676 ztrans_new = 0._wp 667 DO jkb = 1, ipkb ! calculate transport on input grid668 ztrans = ztrans + pdta_read(jb,1,jkb) * pdta_read_dz(jb, 1,jkb)669 ENDDO670 677 DO jk = 1, jpk ! calculate transport on model grid 671 678 ztrans_new = ztrans_new + pdta(jb,1,jk ) * e3u(ji,jj,jk,Kmm ) * umask(ji,jj,jk) … … 674 681 IF(ldtotvel) THEN ! bdy data are total velocity so adjust bt transport term to match input data 675 682 pdta(jb,1,jk) = ( pdta(jb,1,jk) + ( ztrans - ztrans_new ) * r1_hu(ji,jj,Kmm) ) * umask(ji,jj,jk) 676 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero683 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 677 684 pdta(jb,1,jk) = pdta(jb,1,jk) + ( 0._wp - ztrans_new ) * r1_hu(ji,jj,Kmm) * umask(ji,jj,jk) 678 685 ENDIF 679 686 ENDDO 680 687 ENDDO 681 ENDIF 682 683 IF(kgrd == 3) THEN ! do we need to adjust the transport term? 688 689 CASE(3) ! adjust the V-transport term 684 690 DO jb = 1, ipi 685 691 ji = idx_bdy(kbdy)%nbi(jb,kgrd) 686 692 jj = idx_bdy(kbdy)%nbj(jb,kgrd) 687 zh = SUM(pdta_read_dz(jb,1,:) )688 693 ztrans = 0._wp 694 DO jkb = 1, ipkb ! calculate transport on input grid 695 IF( pdta_read(jb,1,jkb) /= pfv ) ztrans = ztrans + pdta_read(jb,1,jkb) * pdta_read_dz(jb,1,jkb) 696 ENDDO 689 697 ztrans_new = 0._wp 690 DO jkb = 1, ipkb ! calculate transport on input grid691 ztrans = ztrans + pdta_read(jb,1,jkb) * pdta_read_dz(jb, 1,jkb)692 ENDDO693 698 DO jk = 1, jpk ! calculate transport on model grid 694 ztrans_new = ztrans_new + pdta(jb,1,jk ) * 699 ztrans_new = ztrans_new + pdta(jb,1,jk ) * e3v(ji,jj,jk,Kmm ) * vmask(ji,jj,jk) 695 700 ENDDO 696 701 DO jk = 1, jpk ! make transport correction 697 702 IF(ldtotvel) THEN ! bdy data are total velocity so adjust bt transport term to match input data 698 703 pdta(jb,1,jk) = ( pdta(jb,1,jk) + ( ztrans - ztrans_new ) * r1_hv(ji,jj,Kmm) ) * vmask(ji,jj,jk) 699 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero704 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 700 705 pdta(jb,1,jk) = pdta(jb,1,jk) + ( 0._wp - ztrans_new ) * r1_hv(ji,jj,Kmm) * vmask(ji,jj,jk) 701 706 ENDIF 702 707 ENDDO 703 708 ENDDO 704 END IF709 END SELECT 705 710 706 711 END SUBROUTINE fld_bdy_interp
Note: See TracChangeset
for help on using the changeset viewer.
