New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 5891 for branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP – NEMO

Ignore:
Timestamp:
2015-11-17T11:15:01+01:00 (8 years ago)
Author:
jamesharle
Message:

bugfix in bdy_interp

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r5705 r5891  
    789789 
    790790      IF( PRESENT(jpk_bdy) .AND. jpk_bdy>0 ) THEN       ! boundary data not on model grid: vertical interpolation 
     791         CALL iom_getatt(num, '_FillValue', fv, cdvar=clvar ) 
     792         dta_read(:,:,:) = -ABS(fv) 
     793         dta_read_z(:,:,:) = 0._wp 
     794         dta_read_dz(:,:,:) = 0._wp 
    791795         CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:jpk_bdy), nrec ) 
    792796         SELECT CASE( igrd )                   
    … …  
    801805               CALL iom_get ( num, jpdom_unknown, 'e3v', dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) ) 
    802806         END SELECT 
    803          CALL iom_getatt(num, '_FillValue', fv, cdvar=clvar ) 
    804807 
    805808#if defined key_bdy 
    … …  
    868871       
    869872      fv_alt = -ABS(fv)  ! set _FillValue < 0 as we make use of MAXVAL and MAXLOC later 
    870        
     873      DO ib = 1, ipi 
     874            zij = idx_bdy(ibdy)%nbi(ib,igrd) 
     875            zjj = idx_bdy(ibdy)%nbj(ib,igrd) 
     876         IF(narea==2) WRITE(*,*) 'MAPI', ib, igrd, map%ptr(ib), narea-1, zij, zjj 
     877      ENDDO 
    871878      ! 
    872879      IF ( map%ll_unstruc ) THEN                            ! unstructured open boundary data file 
    … …  
    885892            zjj = idx_bdy(ibdy)%nbj(ib,igrd) 
    886893            zh  = SUM(dta_read_dz(map%ptr(ib),1,:) ) 
     894            ! Warnings to flag differences in the input and model topgraphy - is this useful/necessary? 
     895            SELECT CASE( igrd )                          
     896               CASE(1) 
     897                  IF( ABS( (zh - ht_0(zij,zjj)) / ht_0(zij,zjj)) * tmask(zij,zjj,1) > 0.01_wp ) THEN 
     898                     WRITE(ibstr,"(I10.10)") map%ptr(ib)  
     899                     CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 
     900                 !   IF(lwp) WRITE(*,*) 'DEPTHT', zh, sum(e3t_0(zij,zjj,:), mask=tmask(zij,zjj,:)==1),  ht_0(zij,zjj), map%ptr(ib), ib, zij, zjj 
     901                  ENDIF 
     902               CASE(2) 
     903                  IF( ABS( (zh - hu_0(zij,zjj)) / hu_0(zij,zjj)) * umask(zij,zjj,1) > 0.01_wp ) THEN 
     904                     WRITE(ibstr,"(I10.10)") map%ptr(ib)  
     905                     CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 
     906                     IF(lwp) WRITE(*,*) 'DEPTHU', zh, sum(e3u_0(zij,zjj,:), mask=umask(zij,zjj,:)==1),  sum(umask(zij,zjj,:)), & 
     907                       &                hu_0(zij,zjj), map%ptr(ib), ib, zij, zjj, narea-1  , & 
     908                        &                dta_read(map%ptr(ib),1,:) 
     909                  ENDIF 
     910               CASE(3) 
     911                  IF( ABS( (zh - hv_0(zij,zjj)) / hv_0(zij,zjj)) * vmask(zij,zjj,1) > 0.01_wp ) THEN 
     912                     WRITE(ibstr,"(I10.10)") map%ptr(ib)  
     913                     CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 
     914                  ENDIF 
     915            END SELECT 
    887916            DO ik = 1, ipk                       
    888917               SELECT CASE( igrd )                        
    … …  
    9881017            zjj = idx_bdy(ibdy)%nbj(ib,igrd) 
    9891018            zh  = SUM(dta_read_dz(ji,jj,:) ) 
     1019            ! Warnings to flag differences in the input and model topgraphy - is this useful/necessary? 
     1020            SELECT CASE( igrd )                          
     1021               CASE(1) 
     1022                  IF( ABS( (zh - ht_0(zij,zjj)) / ht_0(zij,zjj)) * tmask(zij,zjj,1) > 0.01_wp ) THEN 
     1023                     WRITE(ibstr,"(I10.10)") map%ptr(ib)  
     1024                     CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 
     1025                 !   IF(lwp) WRITE(*,*) 'DEPTHT', zh, sum(e3t_0(zij,zjj,:), mask=tmask(zij,zjj,:)==1),  ht_0(zij,zjj), map%ptr(ib), ib, zij, zjj 
     1026                  ENDIF 
     1027               CASE(2) 
     1028                  IF( ABS( (zh - hu_0(zij,zjj)) / hu_0(zij,zjj)) * umask(zij,zjj,1) > 0.01_wp ) THEN 
     1029                     WRITE(ibstr,"(I10.10)") map%ptr(ib)  
     1030                     CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 
     1031                  ENDIF 
     1032               CASE(3) 
     1033                  IF( ABS( (zh - hv_0(zij,zjj)) / hv_0(zij,zjj)) * vmask(zij,zjj,1) > 0.01_wp ) THEN 
     1034                     WRITE(ibstr,"(I10.10)") map%ptr(ib)  
     1035                     CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 
     1036                  ENDIF 
     1037            END SELECT 
    9901038            DO ik = 1, ipk                       
    9911039               SELECT CASE( igrd )                          ! coded for sco - need zco and zps option using min 
    … …  
    10241072 
    10251073         IF(igrd == 2) THEN                                 ! do we need to adjust the transport term? 
    1026             jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 
    1027             ji=map%ptr(ib)-(jj-1)*ilendta 
    1028             zij = idx_bdy(ibdy)%nbi(ib,igrd) 
    1029             zjj = idx_bdy(ibdy)%nbj(ib,igrd) 
    10301074            DO ib = 1, ipi 
    1031               zh = SUM(dta_read_dz(ji,jj,:) ) 
    1032               ztrans = 0._wp 
    1033               ztrans_new = 0._wp 
    1034               DO ik = 1, jpk_bdy                            ! calculate transport on input grid 
     1075               jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 
     1076               ji=map%ptr(ib)-(jj-1)*ilendta 
     1077               zij = idx_bdy(ibdy)%nbi(ib,igrd) 
     1078               zjj = idx_bdy(ibdy)%nbj(ib,igrd) 
     1079               zh = SUM(dta_read_dz(ji,jj,:) ) 
     1080               ztrans = 0._wp 
     1081               ztrans_new = 0._wp 
     1082               DO ik = 1, jpk_bdy                            ! calculate transport on input grid 
    10351083                  ztrans = ztrans + dta_read(ji,jj,ik) * dta_read_dz(ji,jj,ik) 
    1036               ENDDO 
    1037               DO ik = 1, ipk                                ! calculate transport on model grid 
     1084               ENDDO 
     1085               DO ik = 1, ipk                                ! calculate transport on model grid 
    10381086                  ztrans_new = ztrans_new + dta(ib,1,ik) * e3u_0(zij,zjj,ik) * umask(zij,zjj,ik) 
    1039               ENDDO 
    1040               DO ik = 1, ipk                                ! make transport correction 
    1041                  IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 
    1042                     dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * hur(zij,zjj) ) * umask(zij,zjj,ik) 
    1043                  ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 
    1044                     dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp  - ztrans_new ) * hur(zij,zjj) ) * umask(zij,zjj,ik) 
    1045                  ENDIF 
    1046               ENDDO 
     1087               ENDDO 
     1088               DO ik = 1, ipk                                ! make transport correction 
     1089                  IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 
     1090                     dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * hur(zij,zjj) ) * umask(zij,zjj,ik) 
     1091                  ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 
     1092                     dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp  - ztrans_new ) * hur(zij,zjj) ) * umask(zij,zjj,ik) 
     1093                  ENDIF 
     1094               ENDDO 
    10471095            ENDDO 
    10481096         ENDIF 
    10491097 
    10501098         IF(igrd == 3) THEN                                 ! do we need to adjust the transport term? 
    1051            DO ib = 1, ipi 
    1052               jj  = 1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 
    1053               ji  = map%ptr(ib)-(jj-1)*ilendta 
    1054               zij = idx_bdy(ibdy)%nbi(ib,igrd) 
    1055               zjj = idx_bdy(ibdy)%nbj(ib,igrd) 
    1056               zh  = SUM(dta_read_dz(ji,jj,:) ) 
    1057               ztrans = 0._wp 
    1058               ztrans_new = 0._wp 
    1059               DO ik = 1, jpk_bdy                            ! calculate transport on input grid 
     1099            DO ib = 1, ipi 
     1100               jj  = 1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 
     1101               ji  = map%ptr(ib)-(jj-1)*ilendta 
     1102               zij = idx_bdy(ibdy)%nbi(ib,igrd) 
     1103               zjj = idx_bdy(ibdy)%nbj(ib,igrd) 
     1104               zh  = SUM(dta_read_dz(ji,jj,:) ) 
     1105               ztrans = 0._wp 
     1106               ztrans_new = 0._wp 
     1107               DO ik = 1, jpk_bdy                            ! calculate transport on input grid 
    10601108                  ztrans     = ztrans     + dta_read(ji,jj,ik) * dta_read_dz(ji,jj,ik) 
    1061               ENDDO 
    1062               DO ik = 1, ipk                                ! calculate transport on model grid 
     1109               ENDDO 
     1110               DO ik = 1, ipk                                ! calculate transport on model grid 
    10631111                  ztrans_new = ztrans_new + dta(ib,1,ik) * e3v_0(zij,zjj,ik) * vmask(zij,zjj,ik) 
    1064               ENDDO 
    1065               DO ik = 1, ipk                                ! make transport correction 
    1066                  IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 
    1067                     dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * hvr(zij,zjj) ) * vmask(zij,zjj,ik) 
    1068                  ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 
    1069                     dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp  - ztrans_new ) * hvr(zij,zjj) ) * vmask(zij,zjj,ik) 
    1070                  ENDIF 
    1071               ENDDO 
     1112               ENDDO 
     1113               DO ik = 1, ipk                                ! make transport correction 
     1114                  IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 
     1115                     dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * hvr(zij,zjj) ) * vmask(zij,zjj,ik) 
     1116                  ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 
     1117                     dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp  - ztrans_new ) * hvr(zij,zjj) ) * vmask(zij,zjj,ik) 
     1118                  ENDIF 
     1119               ENDDO 
    10721120            ENDDO 
    10731121         ENDIF 
Note: See TracChangeset for help on using the changeset viewer.