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 = 163.0, ywidth = 40.0, 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 == 'rtwn' .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)
  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   ! Pacific removal
    DO j = 1,85
      DO i = 1,69
        fsm(i,j) = 0.
        dum(i,j) = 0.
        dvm(i,j) = 0.
        h(i,j) = 1.
      END DO
    END DO
    DO j = 1,55
      DO i = 70,142
        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