PROGRAM rtofprep

!    Name: rtofprep
!    Subroutines: readrtfinp
!                 outncf
!
!    Author: Biju Thomas on Jan 12, 2016 (GSO-URI)
!
!  1. Purpose :
!
! Read GLOBAL RTOF data (binary) and inerpolate to MPIPOM grids 
! Pepare initial conditions for MPIPOM-TC.
!


USE varsmd
USE outncfmd
USE setupmd
USE fillupmd
USE smooth2d_md
USE readrtfinpmd
IMPLICIT NONE

  INTEGER :: ionedim, allo_err
  INTEGER :: i, j, k, nflg, igfssst
  INTEGER, PARAMETER :: maskval = 0
  INTEGER, DIMENSION(:,:), ALLOCATABLE :: mask
  REAL :: blon, blat, whm
  REAL, PARAMETER :: d2r = 3.1416/180., rearth = 6371.e3
  REAL, PARAMETER :: xwidth = 83.2, ywidth = 37.5, hmax = 5500.0
  REAL, PARAMETER :: largetval = 1.2676506E+30, largewval=1.2676506E+30
  REAL, DIMENSION(:,:), ALLOCATABLE :: f2d
  CHARACTER(LEN=4) :: odata

  READ(5,*) odata ! Region-specific information
  WRITE(6,*) odata
  CALL swcorner(odata, blon, blat)
  WRITE(6,*) blon, blat
  CALL modeldf()

! MPIPOM-TC full (z) and half (zz) sigma levels
  DO k = 1,kb
    z(k) = -z1(k)/z1(kb)
  END DO
  DO k = 1,kb-1
    zz(k) = 0.5*(z(k)+z(k+1))
  END DO
  zz(kb) = z(kb)

  DO j = 1,jm
    DO i = 1,im
      east_e(i,j) = blon+(i-1)*xwidth/(im-1)
      north_e(i,j) = blat+(j-1)*ywidth/(jm-1)
    END DO
  END DO

! MPIPOM-TC grid spacing

  DO j = 1,jm
    DO i = 2,im-1
      dx(i,j) = 0.5*d2r*rearth*sqrt(((east_e(i+1,j)-east_e(i-1,j))* &
        cos(north_e(i,j)*d2r))**2.+(north_e(i+1,j)-north_e(i-1,j))**2.)
    END DO
    dx(1,j) = dx(2,j)
    dx(im,j) = dx(im-1,j)
  END DO

  DO i = 1,im
    DO j = 2,jm-1
      dy(i,j) = 0.5*d2r*rearth*sqrt(((east_e(i,j+1)-east_e(i,j-1))* &
        cos(north_e(i,j)*d2r))**2.+(north_e(i,j+1)-north_e(i,j-1))**2.)
    END DO
    dy(i,1) = dy(i,2)
    dy(i,jm) = dy(i,jm-1)
  END DO

! MPIPOM-TC cell corners

  DO j = 2,jm
    DO i = 2,im
      east_c(i,j) = 0.25*(east_e(i,j)+east_e(i,j-1)+ &
                         east_e(i-1,j)+east_e(i-1,j-1))
      north_c(i,j) = 0.25*(north_e(i,j)+north_e(i,j-1)+ &
                          north_e(i-1,j)+north_e(i-1,j-1))
    END DO
  END DO
  DO i = 2,im
    east_c(i,1) = 2.*east_c(i,2)-east_c(i,3)
    north_c(i,1) = 2.*north_c(i,2)-north_c(i,3)
  END DO
  east_c(1,1) = 2.*east_c(2,1)-east_c(3,1)
  DO j = 2,jm
    east_c(1,j) = 2.*east_c(2,j)-east_c(3,j)
    north_c(1,j) = 2.*north_c(2,j)-north_c(3,j)
  END DO
  north_c(1,1) = 2.*north_c(1,2)-north_c(1,3)

  READ(5,*) ionedim
  IF (ionedim.eq.1) THEN

! MPIPOM-TC u points

    DO i = 1,im
      DO j = 1,jm
        east_u(i,j) = east_e(i,j)
        north_u(i,j) = north_e(i,j)
      END DO
    END DO

! MPIPOM-TC v points

    DO j = 1,jm
      DO i = 1,im
        east_v(i,j) = east_e(i,j)
        north_v(i,j) = north_e(i,j)
      END DO
    END DO

  ELSE

! MPIPOM-TC u points

    DO i = 1,im
      DO j = 1,jm-1
        east_u(i,j) = 0.5*(east_c(i,j)+east_c(i,j+1))
        north_u(i,j) = 0.5*(north_c(i,j)+north_c(i,j+1))
       END DO
      east_u(i,jm) = 0.5*(east_c(i,jm)*3.-east_c(i,jm-1))
      north_u(i,jm) = 0.5*(north_c(i,jm)*3.-north_c(i,jm-1))
    END DO

! MPIPOM-TC v points

    DO j = 1,jm
      DO i = 1,im-1
        east_v(i,j) = 0.5*(east_c(i,j)+east_c(i+1,j))
        north_v(i,j) = 0.5*(north_c(i,j)+north_c(i+1,j))
      END DO
      east_v(im,j) = 0.5*(east_c(im,j)*3.-east_c(im-1,j))
      north_v(im,j) = 0.5*(north_c(im,j)*3.-north_c(im-1,j))
    END DO

  END IF

! MPIPOM-TC rotation angle

  DO j = 1,jm
    DO i = 1,im-1
      rot(i,j) = atan2(north_e(i+1,j)-north_e(i,j), &
                       east_e(i+1,j)-east_e(i,j))
    END DO
    rot(im,j) = rot(im-1,j)
  END DO

  IF (blon < 0.0 .AND. (odata == 'rttr' .OR. odata == 'rtep' .OR. &
      odata == 'rtse' .OR. odata == 'rtsa') ) THEN
     whm = 360.0
  ELSE
     whm = 0.0
  ENDIF 
  CALL readrtfinp(im, jm, nl, east_e+whm, north_e,east_u+whm, north_u,     &
                  east_v+whm, north_v, zrtof, t, s, uin, vin, h)
  write(15)h
  write(8,*)h(1003,28)
  ALLOCATE(mask(im,jm), STAT=allo_err)
  DO j = 1,jm
    DO i = 1,im
      IF ( h(i,j) == largetval ) THEN
         mask(i,j) = maskval
      ELSE
         mask(i,j) = 1
      ENDIF
    ENDDO
  ENDDO
  write(8,*)'a',h(1003,28)
  CALL  smooth2d9ptwave(maskval,mask, h, im, jm)
  write(8,*)'b',h(1003,28)
  CALL  smooth2d9ptwave(maskval,mask, h, im, jm)
  write(8,*)'c',h(1003,28) 
  CALL  smooth2d9ptwave(maskval,mask, h, im, jm)
  write(8,*)'d',h(1003,28)
  IF ( odata == 'rtcp' ) THEN
    CALL  smooth2d9ptwave(maskval,mask, h, im, jm)
    CALL  smooth2d9ptwave(maskval,mask, h, im, jm)
    CALL  smooth2d9ptwave(maskval,mask, h, im, jm)
    CALL  smooth2d9ptwave(maskval,mask, h, im, jm)
    CALL  smooth2d9ptwave(maskval,mask, h, im, jm)
    CALL  smooth2d9ptwave(maskval,mask, h, im, jm)
    CALL  smooth2d9ptwave(maskval,mask, h, im, jm)
    CALL  smooth2d9ptwave(maskval,mask, h, im, jm)
    CALL  smooth2d9ptwave(maskval,mask, h, im, jm)
  END IF

  DEALLOCATE(mask)
  write(8,*)h(1003,28)

!  IF ( .NOT. ALLOCATED(f2d)) ALLOCATE(f2d(im,jm))
!  f2d = h
!  DO j = 1, jm
!    DO i = 1, im
!      DO k = 2, nl
!        IF( zrtof(k) <=  h(i,j) .AND. h(i,j) /= largetval) THEN
!          IF ( t(i,j,k) == largetval .OR. s(i,j,k) == largetval ) THEN
!            f2d(i,j) = zrtof(k-1)
!            EXIT
!          ENDIF
!        ENDIF
!      ENDDO
!    ENDDO
!  ENDDO
!  h = f2d
!  DEALLOCATE(f2d)

  CALL fillup(largetval, zrtof, h, t, s)

  DO j = 1, jm
    DO i = 1, im
       DO k = 1, nl
         IF( zrtof(k) >  h(i,j) ) THEN
            t(i,j,k) = largetval
            s(i,j,k) = largetval
            uin(i,j,k) = largetval
            vin(i,j,k) = largetval
         ENDIF
       ENDDO
    ENDDO
  ENDDO

! Set uin, vin, and el to 0

  DO k = 1,nl
    DO j = 1,jm
      DO i = 1,im
        uin(i,j,k) = 0. ! calculated as geostrophic in old POM-TC
        vin(i,j,k) = 0. ! calculated as geostrophic in old POM-TC
      END DO
    END DO
  END DO

  DO j = 1,jm
    DO i = 1,im
      el(i,j) = 0.
    END DO
  END DO
 
! Neutral profile generation
  DO j = 1,jm
    DO i = 1,im
      nflg = 0
      DO k = 2,nl
        IF (t(i,j,k) == largetval .OR. s(i,j,k) ==  largetval .OR. nflg .EQ. 1) THEN
          t(i,j,k) = t(i,j,k-1)
          s(i,j,k) = s(i,j,k-1)
          nflg = 1
        END IF
      END DO
    END DO
  END DO


! MPIPOM-TC land/sea mask

  write(8,*)h(1003,28)
  DO j = 1,jm
    DO i = 1,im
      IF (h(i,j) < 11.) h(i,j) = 1.  ! Restrict shallowest topography to 11-m depth
      IF (h(i,j) > 1. .AND. h(i,j) /= largetval) THEN
        h(i,j) = MIN(h(i,j), hmax)  ! Restrict bottom topography to 5500-m depth
        fsm(i,j) = 1.
        dum(i,j) = 1.
        dvm(i,j) = 1.
      ELSE
        fsm(i,j) = 0.
        dum(i,j) = 0.
        dvm(i,j) = 0.
        h(i,j) = 1.
      END IF
    END DO
  END DO


  IF (odata == 'rttr') THEN    

    DO j = 1, 110
    DO i = 1, 100
       fsm(i,j) = 0.
        dum(i,j) = 0.
        dvm(i,j) = 0.
        h(i,j) = 1.
    END DO
    END DO

    DO j = 1, 75
    DO i = 100, 150
       fsm(i,j) = 0.
        dum(i,j) = 0.
        dvm(i,j) = 0.
        h(i,j) = 1.
    END DO
    END DO

    DO j = 1, 25
    DO i = 150, 167
       fsm(i,j) = 0.
        dum(i,j) = 0.
        dvm(i,j) = 0.
        h(i,j) = 1.
    END DO
    END DO

    DO j = 1, 16
    DO i = 167, 196
       fsm(i,j) = 0.
        dum(i,j) = 0.
        dvm(i,j) = 0.
        h(i,j) = 1.
    END DO
    END DO

    DO j = 1, 22
    DO i = 196, 215
       fsm(i,j) = 0.
        dum(i,j) = 0.
        dvm(i,j) = 0.
        h(i,j) = 1.
    END DO
    END DO

    DO j = 1, 2
    DO i = 215, 235
       fsm(i,j) = 0.
        dum(i,j) = 0.
        dvm(i,j) = 0.
        h(i,j) = 1.
    END DO
    END DO
  END IF

  IF (odata == 'rtep') THEN   ! Caribbean and Gulf of America removal
    DO j = 120,jm
    DO i = 810,im
      fsm(i,j) = 0.
      dum(i,j) = 0.
      dvm(i,j) = 0.
      h(i,j) = 1.
    END DO
    END DO

    DO j = 145,jm
    DO i = 715,im
      fsm(i,j) = 0.
      dum(i,j) = 0.
      dvm(i,j) = 0.
      h(i,j) = 1.
    END DO
    END DO
  END IF

  IF (odata == 'rtsa') THEN   ! Southeast Pacific removal
    DO j = 1,jm
      DO i = 1,194
        fsm(i,j) = 0.
        dum(i,j) = 0.
        dvm(i,j) = 0.
        h(i,j) = 1.
      END DO
    END DO
  END IF

  IF (ionedim.ne.1) THEN
! MPIPOM-TC land/sea mask modification
    DO j = 1,jm-1
      DO i = 1,im
        IF (fsm(i,j) == 0. .AND. fsm(i,j+1) /= 0.) dvm(i,j+1) = 0.
      END DO
    END DO

    DO j = 1,jm
      DO i = 1,im-1
        IF (fsm(i,j) == 0. .AND. fsm(i+1,j) /= 0.) dum(i+1,j) = 0. 
      END DO
    END DO
  END IF

! Read and assimilate gfs sst into temperature field

 READ(5,*) igfssst
 WRITE(6,*) 'SST assimilation (0=no; 1=yes)? ',igfssst
! IF (igfssst.eq.1) THEN
!     print*, 'Preparing to read gfs sst...'
!     read(23,224) im1,jm1 ! ensure lonlat.gfs   is linked to fort.23
!     rewind(23)           ! ensure sst.gfs.dat  is linked to fort.21
! 224 format(1x,2i5)       ! ensure mask.gfs.dat is linked to fort.22
!     call serftempr(im1,jm1,stm,fsm,alon,alat,im,jm)
!     print*, '... serftempr read gfs sst: im1,jm = ',im1,jm1
!
!     print*, 'Preparing to assimilate gfs sst...'
!     call mixsstz(t,stm,h,fsm,zrtof,im,jm,nl,0)
!     print*, '... mixsstz assimilated gfs sst'
! END IF

! MPIPOM-TC temperature and salinity modification

  DO k = 1,nl
    DO j = 1,jm
       DO i = 1,im
         t(i,j,k) = t(i,j,k)*fsm(i,j)
         s(i,j,k) = s(i,j,k)*fsm(i,j)
         IF( zrtof(k) <=  h(i,j) .AND. fsm(i,j) == 1.0 ) THEN
           IF (t(i,j,k) == largetval .OR. s(i,j,k) == largetval) THEN
             t(i,j,k) = 0.
             s(i,j,k) = 0.
             uin(i,j,k) = 0.
             vin(i,j,k) = 0.
             PRINT*, 'Bprobl',  east_e(i,j),north_e(i,j),k, h(i,j),fsm(i,j)
          END IF
        END IF
      END DO
    END DO
  END DO

! (for zrtof > h replace largetval with neutral profile)

!  DO j = 1,jm
!    DO i = 1,im
!      DO k = 2,nl
!        IF ( zrtof(k) > h(i,j) .AND. fsm(i,j) == 1.0 ) THEN
!          IF (t(i,j,k) == largetval .OR. s(i,j,k) ==  largetval ) THEN
!            t(i,j,k) = t(i,j,k-1)
!            s(i,j,k) = s(i,j,k-1)
!          END IF
!        END IF
!      END DO
!    END DO
!  END DO


! MPIPOM-TC u and v modification

  DO k = 1,nl
    IF (ionedim.eq.1) THEN
      DO j = 1,jm
        DO i = 1,im
          u(i,j,k) = uin(i,j,k)*dum(i,j)
          v(i,j,k) = vin(i,j,k)*dvm(i,j)
        END DO
      END DO
    ELSE
      DO j = 1,jm
        u(1,j,k) = uin(1,j,k)*dum(1,j)
          DO i = 2,im
             u(i,j,k) = 0.5*(uin(i-1,j,k)+uin(i,j,k))*dum(i,j)
          END DO
      END DO
      DO i = 1,im
        v(i,1,k) = vin(i,1,k)*dvm(i,1)
        DO j = 2,jm
          v(i,j,k) = 0.5*(vin(i,j-1,k)+vin(i,j,k))*dvm(i,j)
        END DO
      END DO
    END IF
  END DO

! MPIPOM-TC el modification

 DO j = 1,jm
   DO i = 1,im
     el(i,j) = el(i,j)*fsm(i,j)
   END DO
 END DO

! MPIPOM-TC tclim and sclim calcuation

   DO k = 1,nl
     taa(k) = 0.
     saa(k) = 0.
     area(k) = 0.
     DO j = 1,jm
       DO i = 1,im
         IF (t(i,j,k) /= 0.) THEN 
           taa(k) = taa(k)+t(i,j,k)
           saa(k) = saa(k)+s(i,j,k)
           area(k) = area(k)+1.
         END IF
       END DO
     END DO
     taa(k) = taa(k)/area(k)
     saa(k) = saa(k)/area(k)
   END DO

   DO k = 1,nl
     DO j = 1,jm
       DO i = 1,im
          tclim(i,j,k) = taa(k)
          sclim(i,j,k) = saa(k)
       END DO
     END DO
   END DO

! Prepare to create MPIPOM-TC input files

  READ(5,*) runid
  WRITE(6,*) TRIM(ADJUSTL(runid))
  write(8,*)h(1003,28)
  
  CALL outncf

END PROGRAM