Changeset 5625

Timestamp:

2015-07-22T09:48:42+02:00 (9 years ago)

Author:

jamesharle

Message:

Update interpolation code to handle structured BDY.

File:

• branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 (modified) (10 diffs)

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

@@ -773 +773 @@
          SELECT CASE( igrd )                  
             CASE(1)
+               WRITE(*,*) 'reading gdept' 
                CALL iom_get ( num, jpdom_unknown, 'gdept', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) )
+               WRITE(*,*) 'reading e3t' 
                CALL iom_get ( num, jpdom_unknown, 'e3t', dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) )
+               WRITE(*,*) 'finished reading' 
             CASE(2)  
                CALL iom_get ( num, jpdom_unknown, 'gdepu', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) )
@@ -785 +788 @@
 
 #if defined key_bdy
+      WRITE(*,*) 'going in bdy interp' 
          CALL fld_bdy_interp(dta_read, dta_read_z, dta_read_dz, map, jpk_bdy, igrd, ibdy, fv, dta, fvl)
+      WRITE(*,*) 'coming out bdy interp' 
 #endif
       ELSE ! boundary data assumed to be on model grid
@@ -839 +844 @@
       REAL(wp)                                ::   fv_alt, ztrans, ztrans_new ! fillvalue and alternative -ABS(fv)
       !!---------------------------------------------------------------------
+     
+      WRITE(*,*) 'in bdy interp' 
 
       ipi       = SIZE( dta, 1 )
@@ -847 +854 @@
  
       fv_alt = -ABS(fv)  ! set _FillValue < 0 as we make use of MAXVAL and MAXLOC later
+      
       !
       IF ( map%ll_unstruc ) THEN ! unstructured open boundary data file
+
          DO ib = 1, ipi
             DO ik = 1, jpk_bdy
@@ -879 +888 @@
 
          IF(igrd == 2) THEN ! do we need to adjust the transport term?
-           DO ib = 1, ipi
+            DO ib = 1, ipi
               zh = SUM(dta_read_dz(map%ptr(ib),1,:) )
               ztrans = 0._wp
@@ -892 +901 @@
               DO ik = 1, ipk
                  zdz =  e3u_0(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd),ik)   
-                 zz  =  hur(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd))   
+                 zz  =  hur(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd))       ! Is array HUR the correct one to use????
                  IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data
                     dta(ib,1,ik) = dta(ib,1,ik) + ( ztrans - ztrans_new ) * ( zdz * zz )
@@ -903 +912 @@
 
          IF(igrd == 3) THEN ! do we need to adjust the transport term?
-           DO ib = 1, ipi
+            DO ib = 1, ipi
               zh = SUM(dta_read_dz(map%ptr(ib),1,:) )
               ztrans = 0._wp
@@ -931 +940 @@
 
       ELSE ! structured open boundary file
+
          DO ib = 1, ipi
             jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta))
@@ -939 +949 @@
                   dta_read_dz(ji,jj,ik) = 0._wp ! safety: put 0._wp into external thickness factors to ensure transport is correct
                ENDIF
-     !            dta(ib,1,ik) = fv_alt    ! put fillvalue into new field as if all goes well all wet points will be replaced
             ENDDO
          ENDDO 
-      !
+       
          DO ib = 1, ipi
             jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta))
@@ -964 +973 @@
             END DO
          END DO
-      ENDIF
+
+         IF(igrd == 2) THEN ! do we need to adjust the transport term?
+            jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta))
+            ji=map%ptr(ib)-(jj-1)*ilendta
+            DO ib = 1, ipi
+              zh = SUM(dta_read_dz(ji,jj,:) )
+              ztrans = 0._wp
+              ztrans_new = 0._wp
+              DO ik = 1, jpk_bdy
+                  ztrans = ztrans + dta_read(ji,jj,ik) * dta_read_dz(ji,jj,ik)
+              ENDDO
+              DO ik = 1, ipk
+                  zdz =  e3u_0(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd),ik)
+                  ztrans_new = ztrans_new + dta(ib,1,ik) * zdz
+              ENDDO
+              DO ik = 1, ipk
+                 zdz =  e3u_0(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd),ik)
+                 zz  =  hur(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd))
+                 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data
+                    dta(ib,1,ik) = dta(ib,1,ik) + ( ztrans - ztrans_new ) * ( zdz * zz )
+                 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero
+                    dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * ( zdz * zz )
+                 ENDIF
+              ENDDO
+            ENDDO
+         ENDIF
+
+         IF(igrd == 3) THEN ! do we need to adjust the transport term?
+            jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta))
+            ji=map%ptr(ib)-(jj-1)*ilendta
+           DO ib = 1, ipi
+              zh = SUM(dta_read_dz(ji,jj,:) )
+              ztrans = 0._wp
+              ztrans_new = 0._wp
+              DO ik = 1, jpk_bdy
+                  ztrans = ztrans + dta_read(ji,jj,ik) * dta_read_dz(ji,jj,ik)
+              ENDDO
+              DO ik = 1, ipk
+                  zdz =  e3v_0(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd),ik)
+                  ztrans_new = ztrans_new + dta(ib,1,ik) * zdz
+              ENDDO
+              DO ik = 1, ipk
+                 zdz =  e3v_0(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd),ik)
+                 zz  =  hvr(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd))
+                 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data
+                    dta(ib,1,ik) = dta(ib,1,ik) + ( ztrans - ztrans_new ) * ( zdz * zz )
+                 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero
+                    dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * ( zdz * zz )
+                 ENDIF
+              ENDDO
+            ENDDO
+         ENDIF
+
+      ENDIF ! endif unstructured or structured
 
    END SUBROUTINE fld_bdy_interp