source: branches/2015/dev_r5187_UKMO13_simplification/NEMOGCM/TOOLS/SCOORD_GEN/src/utils.F90 @ 5794

MODULE utils

   USE netcdf

   IMPLICIT NONE
   PUBLIC             ! allows the acces to par_oce when dom_oce is used
   !                  ! exception to coding rules... to be suppressed ???

!   PUBLIC dom_oce_alloc

   INTEGER, PARAMETER   :: dp=8 , sp=4, wp=dp

   !! All coordinates
   !! ---------------
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   gdep3w_0           !: depth of t-points (sum of e3w) (m)
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   gdept_0, gdepw_0   !: analytical (time invariant) depth at t-w  points (m)
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   gphit              !: latitude at t points
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e3v_0  , e3f_0     !: analytical (time invariant) vertical scale factors at  v-f
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e3t_0  , e3u_0     !:                                      t-u  points (m)
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e3vw_0             !: analytical (time invariant) vertical scale factors at  vw
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e3w_0  , e3uw_0    !:                                      w-uw points (m)

   !! s-coordinate and hybrid z-s-coordinate
   !! =----------------======---------------
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hbatv , hbatf      !: ocean depth at the vertical of  v--f
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hbatt , hbatu      !:                                 t--u points (m)
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   scosrf, scobot     !: ocean surface and bottom topographies 
   !                                                                           !  (if deviating from coordinate surfaces in HYBRID)
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hifv  , hiff       !: interface depth between stretching at v--f
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hift  , hifu       !: and quasi-uniform spacing             t--u points (m)
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rx1                !: Maximum grid stiffness ratio

   !!----------------------------------------------------------------------
   !! masks, bathymetry
   !! ---------------------------------------------------------------------
   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbathy             !: number of ocean level (=0, 1, ... , jpk-1)
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bathy              !: ocean depth (meters - read from file)

   !! Other variables needed by scoord_gen
   INTEGER  ::   jpi, jpj, jpk            ! Size of the domain - read from bathy or namelist?
   INTEGER  ::   ji, jj, jk, jl           ! dummy loop argument
   INTEGER  ::   iip1, ijp1, iim1, ijm1   ! temporary integers
   INTEGER  ::   ios                      ! Local integer output status for namelist read and allocation
   INTEGER,PARAMETER  ::   numnam=8       ! File handle for namelist 
   REAL(wp) ::   zrmax, ztaper   ! temporary scalars
   REAL(wp) ::   zrfact
   !
   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2
   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zenv, ztmp, zmsk, zri, zrj, zhbat

   !Namelist variables
   REAL(wp) :: rn_jpk, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax, rn_theta
   REAL(wp) :: rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b
   LOGICAL :: ln_s_sh94, ln_s_sf12, ln_sigcrit, ln_eq_taper
   CHARACTER(len=50) :: cn_coord_hgr

   NAMELIST/namzgr_sco/rn_jpk, ln_s_sh94, ln_s_sf12, ln_sigcrit, ln_eq_taper, & 
                  &      cn_coord_hgr, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, &
                  &      rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b

  ! IDs for output netcdf file
  INTEGER :: id_x, id_y, id_z
  INTEGER :: ncout
  INTEGER, DIMENSION(20) :: var_ids  !Array to contain all variable IDs

   CONTAINS

   INTEGER FUNCTION dom_oce_alloc()
      !!----------------------------------------------------------------------
      INTEGER, DIMENSION(4) :: ierr
      !!----------------------------------------------------------------------
      ierr(:) = 0
      !
      ALLOCATE( zenv(jpi,jpj), ztmp(jpi,jpj), zmsk(jpi,jpj), zri(jpi,jpj), zrj(jpi,jpj), &
         &      zhbat(jpi,jpj) , ztmpi1(jpi,jpj), ztmpi2(jpi,jpj), ztmpj1(jpi,jpj), ztmpj2(jpi,jpj), STAT=ierr(1) )
         !
      ALLOCATE( gdep3w_0(jpi,jpj) , e3v_0(jpi,jpj) , e3f_0(jpi,jpj) ,                         &
         &      gdept_0(jpi,jpj) , e3t_0(jpi,jpj) , e3u_0 (jpi,jpj) ,                         &
         &      gdepw_0(jpi,jpj) , e3w_0(jpi,jpj) , e3vw_0(jpi,jpj) ,                         &
         &      gphit(jpi,jpj)   , e3uw_0(jpi,jpj) , STAT=ierr(2) )
         !
         !
      ALLOCATE( hbatv (jpi,jpj) , hbatf (jpi,jpj) ,     &
         &      hbatt (jpi,jpj) , hbatu (jpi,jpj) ,     &
         &      scosrf(jpi,jpj) , scobot(jpi,jpj) ,     &
         &      hifv  (jpi,jpj) , hiff  (jpi,jpj) ,     &
         &      hift  (jpi,jpj) , hifu  (jpi,jpj) , rx1 (jpi,jpj) , STAT=ierr(3) )

      ALLOCATE( mbathy(jpi,jpj) , STAT=ierr(4) )
     !
      dom_oce_alloc = MAXVAL(ierr)
      !
   END FUNCTION dom_oce_alloc
 

   SUBROUTINE read_bathy()
     !! Read bathymetry from input netcdf file
     INTEGER :: var_id, ncin

     CALL check_nf90( nf90_open('bathy.nc', NF90_NOWRITE, ncin), 'Error opening bathy.nc file' )

     ! Find the size of the input bathymetry
     CALL dimlen(ncin, 'lon', jpi)    
     CALL dimlen(ncin, 'lat', jpj)    
     
     ALLOCATE( bathy(jpi, jpj) )
     
     ! Read the bathymetry variable from file
     CALL check_nf90( nf90_inq_varid( ncin, 'Bathymetry', var_id ), 'Cannot get variable ID for bathymetry')
     CALL check_nf90( nf90_get_var( ncin, var_id, bathy, (/ 1,1 /), (/ jpi, jpj /) ) )

     CALL check_nf90( nf90_close(ncin), 'Error closing bathy.nc file' )

   END SUBROUTINE read_bathy

   SUBROUTINE read_gphit()
   !! Read gphit from horizontal coordinate file if required
     INTEGER :: var_id, ncin

     CALL check_nf90( nf90_open(cn_coord_hgr, NF90_NOWRITE, ncin), 'Error opening horizontal coordinate file' )

     ! Read gphit variable from file
     CALL check_nf90( nf90_inq_varid( ncin, 'gphit', var_id ), 'Cannot get variable ID for bathymetry')
     CALL check_nf90( nf90_get_var( ncin, var_id, gphit, (/ 1,1 /), (/ jpi, jpj /) ) )

     CALL check_nf90( nf90_close(ncin), 'Error closing horizontal coordinate file' )

   END SUBROUTINE read_gphit

   SUBROUTINE dimlen( ncid, dimname, len )
     ! Determine the length of dimension dimname
     INTEGER, INTENT(in)          :: ncid
     CHARACTER(LEN=*), INTENT(in) :: dimname
     INTEGER, INTENT(out)         :: len
     ! Local variables
     INTEGER :: id_var, istatus

     id_var = 1
     CALL check_nf90( nf90_inq_dimid(ncid, dimname, id_var), 'Dimension not found in file')
     CALL check_nf90( nf90_inquire_dimension(ncid,id_var,len=len))

   END SUBROUTINE dimlen
  
 
   SUBROUTINE make_coord_file()
     ! Create new coordinates file and define dimensions and variables ready for
     ! writing
     

     !Create the file
     CALL check_nf90( nf90_create('coord_zgr.nc', NF90_CLOBBER, ncout), 'Could not create output file')
     !
     !Define dimensions
     CALL check_nf90( nf90_def_dim(ncout, 'x', jpi, id_x) )
     CALL check_nf90( nf90_def_dim(ncout, 'y', jpj, id_y) )
     CALL check_nf90( nf90_def_dim(ncout, 'z', jpk, id_z) )
     !
     !Define variables - include all varibles that would be put into the mesh
     !mask file
     CALL check_nf90( nf90_def_var(ncout, 'gdept_0', nf90_double, (/id_x, id_y,id_z/), var_ids(1)) )
     CALL check_nf90( nf90_def_var(ncout, 'gdepw_0', nf90_double, (/id_x, id_y,id_z/), var_ids(2)) )
     CALL check_nf90( nf90_def_var(ncout, 'gdep3w_0', nf90_double, (/id_x, id_y,id_z/), var_ids(3)) )
     CALL check_nf90( nf90_def_var(ncout, 'e3f_0', nf90_double, (/id_x, id_y,id_z/), var_ids(4)) )
     CALL check_nf90( nf90_def_var(ncout, 'e3t_0', nf90_double, (/id_x, id_y,id_z/), var_ids(5)) )
     CALL check_nf90( nf90_def_var(ncout, 'e3u_0', nf90_double, (/id_x, id_y,id_z/), var_ids(6)) )
     CALL check_nf90( nf90_def_var(ncout, 'e3v_0', nf90_double, (/id_x, id_y,id_z/), var_ids(7)) )
     CALL check_nf90( nf90_def_var(ncout, 'e3w_0', nf90_double, (/id_x, id_y,id_z/), var_ids(8)) )
     CALL check_nf90( nf90_def_var(ncout, 'e3uw_0', nf90_double, (/id_x, id_y,id_z/), var_ids(9)) )
     CALL check_nf90( nf90_def_var(ncout, 'e3vw_0', nf90_double, (/id_x, id_y,id_z/), var_ids(10)) )
     ! 2D fields
     CALL check_nf90( nf90_def_var(ncout, 'mbathy', nf90_double, (/id_x, id_y/), var_ids(11)) )
     CALL check_nf90( nf90_def_var(ncout, 'hbatt', nf90_double, (/id_x, id_y/), var_ids(12)) )
     CALL check_nf90( nf90_def_var(ncout, 'hbatu', nf90_double, (/id_x, id_y/), var_ids(13)) )
     CALL check_nf90( nf90_def_var(ncout, 'hbatv', nf90_double, (/id_x, id_y/), var_ids(14)) )
     CALL check_nf90( nf90_def_var(ncout, 'hbatf', nf90_double, (/id_x, id_y/), var_ids(15)) )

     
     ! End define mode
     CALL check_nf90( nf90_enddef(ncout) )
     
     WRITE(*,*) 'Opened coord_zgr.nc file and defined variables'

   END SUBROUTINE make_coord_file

   SUBROUTINE write_netcdf_2d_vars()

     CALL check_nf90( nf90_put_var(ncout, var_ids(11), mbathy, (/ 1,1 /), (/ jpi, jpj /) ) )
     CALL check_nf90( nf90_put_var(ncout, var_ids(12), hbatt, (/ 1,1 /), (/ jpi, jpj /) ) )
     CALL check_nf90( nf90_put_var(ncout, var_ids(13), hbatu, (/ 1,1 /), (/ jpi, jpj /) ) )
     CALL check_nf90( nf90_put_var(ncout, var_ids(14), hbatv, (/ 1,1 /), (/ jpi, jpj /) ) )
     CALL check_nf90( nf90_put_var(ncout, var_ids(15), hbatf, (/ 1,1 /), (/ jpi, jpj /) ) )

   END SUBROUTINE write_netcdf_2d_vars

   SUBROUTINE write_netcdf_3d_vars(kk)
   ! Write  variables to the netcdf file at level kk
     INTEGER, INTENT(in) :: kk

     CALL check_nf90( nf90_put_var(ncout, var_ids(1), gdept_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
     CALL check_nf90( nf90_put_var(ncout, var_ids(2), gdepw_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
     CALL check_nf90( nf90_put_var(ncout, var_ids(3), gdep3w_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
     CALL check_nf90( nf90_put_var(ncout, var_ids(4), e3f_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
     CALL check_nf90( nf90_put_var(ncout, var_ids(5), e3t_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
     CALL check_nf90( nf90_put_var(ncout, var_ids(6), e3u_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
     CALL check_nf90( nf90_put_var(ncout, var_ids(7), e3v_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
     CALL check_nf90( nf90_put_var(ncout, var_ids(8), e3w_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
     CALL check_nf90( nf90_put_var(ncout, var_ids(9), e3uw_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
     CALL check_nf90( nf90_put_var(ncout, var_ids(10), e3vw_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )

   END SUBROUTINE write_netcdf_3d_vars

   SUBROUTINE check_nf90( istat, message )
      !Check for netcdf errors
      INTEGER, INTENT(in) :: istat
      CHARACTER(LEN=*), INTENT(in), OPTIONAL :: message

      IF (istat /= nf90_noerr) THEN
         WRITE(*,*) 'ERROR! : '//TRIM(nf90_strerror(istat))
         IF ( PRESENT(message) ) THEN ; WRITE(*,*) message ; ENDIF
         STOP
      ENDIF

   END SUBROUTINE check_nf90


END MODULE utils