! NAME: WRITE_NC_phys
!
! INPUT:
!       group = group code (ex: IPSL)
!       production = production (model and version ex: OPA 8.1)
!
!       secs_per_month = number of seconds in each model month
!
!       imt_trx = tracer Zonal grid dimension
!       jmt_trx = tracer Meridional grid dimension
!       kmt_trx = tracer Depth grid dimension
!
!       PT = monthly potential temperature [C]
!       SAL = monthly salinity [psu]
!       TFLX = monthly net surface heat flux [W/m^2]
!       H2OFLX = monthly net surface freshwater flux [m/y]
!
!       kmt_UV = Horizontal Velocity Depth grid dimension
!
!       imt_Eul_U = Eulerian Zonal Velocity Zonal grid dimension
!       jmt_Eul_U = Eulerian Zonal Velocity Meridional grid dimension
!       Eul_U = monthly Eulerian Zonal Velocity [m/s]
!       WS_x = monthly Zonal Wind Stress [N/m^2]
!
!       imt_Eul_V = Eulerian Meridional Velocity Zonal grid dimension
!       jmt_Eul_V = Eulerian Meridional Velocity Meridional grid dimension
!       Eul_V = monthly Eulerian Meridional Velocity [m/s]
!       WS_y = monthly Meridional Wind Stress [N/m^2]
!
!       has_eddy = logical flag, true if model has eddy induced velocities
!
!       imt_Eddy_U = Eddy Induced Zonal Velocity Zonal grid dimension
!       jmt_Eddy_U = Eddy Induced Zonal Velocity Meridional grid dimension
!       Eddy_U = monthly Eddy Induced Zonal Velocity [m/s]
!
!       imt_Eddy_V = Eddy Induced Meridional Velocity Zonal grid dimension
!       jmt_Eddy_V = Eddy Induced Meridional Velocity Meridional grid dimension
!       Eddy_V = monthly Eddy Induced Meridional Velocity [m/s]
!
! OUTPUT FILE:
!       A netcdf file (clobbed if exists) with
!       the filename group_phys.nc
!
!***********************************************************************

      SUBROUTINE WRITE_NC_phys(
     &   group, production,
     &   secs_per_month,
     &   imt_trx, jmt_trx, kmt_trx,
     &   PT, SAL, TFLX, H2OFLX,
     &   kmt_UV,
     &   imt_Eul_U, jmt_Eul_U, Eul_U, WS_x,
     &   imt_Eul_V, jmt_Eul_V, Eul_V, WS_y,
     &   has_eddy,
     &   imt_Eddy_U, jmt_Eddy_U, Eddy_U,
     &   imt_Eddy_V, jmt_Eddy_V, Eddy_V)

      IMPLICIT NONE
#include "netcdf.inc"

! Arguments
      CHARACTER*(*) group
      CHARACTER*(*) production
      REAL secs_per_month(12)
      INTEGER imt_trx, jmt_trx, kmt_trx

      REAL PT(imt_trx, jmt_trx, kmt_trx, 12)
      REAL SAL(imt_trx, jmt_trx, kmt_trx, 12)
      REAL TFLX(imt_trx, jmt_trx, 12)
      REAL H2OFLX(imt_trx, jmt_trx, 12)

      INTEGER kmt_UV

      INTEGER imt_Eul_U, jmt_Eul_U
      REAL Eul_U(imt_Eul_U, jmt_Eul_U, kmt_UV, 12)
      REAL WS_x(imt_Eul_U, jmt_Eul_U, 12)

      INTEGER imt_Eul_V, jmt_Eul_V
      REAL Eul_V(imt_Eul_V, jmt_Eul_V, kmt_UV, 12)
      REAL WS_y(imt_Eul_V, jmt_Eul_V, 12)

      LOGICAL has_eddy

      INTEGER imt_Eddy_U, jmt_Eddy_U
      REAL Eddy_U(imt_Eddy_U, jmt_Eddy_U, kmt_UV, 12)

      INTEGER imt_Eddy_V, jmt_Eddy_V
      REAL Eddy_V(imt_Eddy_V, jmt_Eddy_V, kmt_UV, 12)

! Constants
      REAL msv
      PARAMETER (msv = -1.E+34)

! Local variables
      CHARACTER*256 filename
      INTEGER NC_ID
      INTEGER STATUS
      INTEGER TWO_DIM
      INTEGER time_DIM
      INTEGER lon_trx_DIM, lat_trx_DIM, depth_trx_DIM
      INTEGER depth_UV_DIM
      INTEGER lon_Eul_U_DIM, lat_Eul_U_DIM
      INTEGER lon_Eul_V_DIM, lat_Eul_V_DIM
      INTEGER lon_Eddy_U_DIM, lat_Eddy_U_DIM
      INTEGER lon_Eddy_V_DIM, lat_Eddy_V_DIM
      INTEGER DIMIDS(4)
      INTEGER time_ID, bounds_time_ID
      INTEGER PT_ID, SAL_ID, TFLX_ID, H2OFLX_ID
      INTEGER Eul_U_ID, WS_x_ID
      INTEGER Eul_V_ID, WS_y_ID
      INTEGER Eddy_U_ID
      INTEGER Eddy_V_ID

      INTEGER I
      REAL TIME(12), bounds_TIME(12,2)

!***********************************************************************
! Build the NetCDF filename
!***********************************************************************
      filename = group//'_phys.nc'

!***********************************************************************
! Open the NetCDF file
!***********************************************************************
      STATUS = NF_CREATE(filename, NF_CLOBBER, NC_ID)
      IF (STATUS .NE. NF_NOERR) THEN
         PRINT *, 'NetCDF error opening ', filename
         PRINT *, NF_STRERROR(STATUS)
         STOP 'Stopped'
      ENDIF

!***********************************************************************
! Define dimensions
!***********************************************************************
      CALL DEF_DIM(NC_ID, 'two', 2, TWO_DIM)

      CALL DEF_DIM(NC_ID, 'time', 12, time_DIM)

      CALL DEF_DIM(NC_ID, 'x_trx', imt_trx, lon_trx_DIM)
      CALL DEF_DIM(NC_ID, 'y_trx', jmt_trx, lat_trx_DIM)
      CALL DEF_DIM(NC_ID, 'z_trx', kmt_trx, depth_trx_DIM)

      CALL DEF_DIM(NC_ID, 'z_UV', kmt_UV, depth_UV_DIM)

      CALL DEF_DIM(NC_ID, 'x_Eul_U', imt_Eul_U, lon_Eul_U_DIM)
      CALL DEF_DIM(NC_ID, 'y_Eul_U', jmt_Eul_U, lat_Eul_U_DIM)

      CALL DEF_DIM(NC_ID, 'x_Eul_V', imt_Eul_V, lon_Eul_V_DIM)
      CALL DEF_DIM(NC_ID, 'y_Eul_V', jmt_Eul_V, lat_Eul_V_DIM)

      IF (has_eddy) THEN
         CALL DEF_DIM(NC_ID, 'x_Eddy_U', imt_Eddy_U, lon_Eddy_U_DIM)
         CALL DEF_DIM(NC_ID, 'y_Eddy_U', jmt_Eddy_U, lat_Eddy_U_DIM)

         CALL DEF_DIM(NC_ID, 'x_Eddy_V', imt_Eddy_V, lon_Eddy_V_DIM)
         CALL DEF_DIM(NC_ID, 'y_Eddy_V', jmt_Eddy_V, lat_Eddy_V_DIM)
      ENDIF ! IF (has_eddy)

!***********************************************************************
! Define variables & attributes
!***********************************************************************

      DIMIDS(1) = time_DIM
      CALL DEF_REAL(NC_ID, 'time', time_ID, 1, DIMIDS)
      CALL ADD_ATT_TEXT(NC_ID, time_ID, 'long_name', 'Time')
      CALL ADD_ATT_TEXT(NC_ID, time_ID, 'units',
     &   'seconds since 0001-01-01 00:00:00')
      CALL ADD_ATT_TEXT(NC_ID, time_ID, 'bounds', 'bounds_time')

      DIMIDS(1) = time_DIM
      DIMIDS(2) = TWO_DIM
      CALL DEF_REAL(NC_ID, 'bounds_time', bounds_time_ID, 2, DIMIDS)

      DIMIDS(1) = lon_trx_DIM
      DIMIDS(2) = lat_trx_DIM
      DIMIDS(3) = depth_trx_DIM
      DIMIDS(4) = time_DIM

      CALL DEF_REAL(NC_ID, 'PT', PT_ID, 4, DIMIDS)
      CALL ADD_ATT_TEXT(NC_ID, PT_ID, 'long_name',
     &   'Potential Temperature')
      CALL ADD_ATT_TEXT(NC_ID, PT_ID, 'units', 'C')
      CALL ADD_ATT_TEXT(NC_ID, PT_ID, 'coordinates',
     &   'lon lat depth time')
      CALL ADD_ATT_REAL(NC_ID, PT_ID, 'missing_value', msv)

      CALL DEF_REAL(NC_ID, 'SAL', SAL_ID, 4, DIMIDS)
      CALL ADD_ATT_TEXT(NC_ID, SAL_ID, 'long_name', 'Salinity')
      CALL ADD_ATT_TEXT(NC_ID, SAL_ID, 'units', 'psu')
      CALL ADD_ATT_TEXT(NC_ID, SAL_ID, 'coordinates',
     &   'lon lat depth time')
      CALL ADD_ATT_REAL(NC_ID, SAL_ID, 'missing_value', msv)

      DIMIDS(1) = lon_trx_DIM
      DIMIDS(2) = lat_trx_DIM
      DIMIDS(3) = time_DIM

      CALL DEF_REAL(NC_ID, 'TFLX', TFLX_ID, 3, DIMIDS)
      CALL ADD_ATT_TEXT(NC_ID, TFLX_ID, 'long_name',
     &   'Net Surface Heat Flux')
      CALL ADD_ATT_TEXT(NC_ID, TFLX_ID, 'units', 'W/m^2')
      CALL ADD_ATT_TEXT(NC_ID, TFLX_ID, 'coordinates',
     &   'lon lat time')
      CALL ADD_ATT_REAL(NC_ID, TFLX_ID, 'missing_value', msv)

      CALL DEF_REAL(NC_ID, 'H2OFLX', H2OFLX_ID, 3, DIMIDS)
      CALL ADD_ATT_TEXT(NC_ID, H2OFLX_ID, 'long_name',
     &   'Net Surface Freshwater Flux')
      CALL ADD_ATT_TEXT(NC_ID, H2OFLX_ID, 'units', 'm/y')
      CALL ADD_ATT_TEXT(NC_ID, H2OFLX_ID, 'coordinates',
     &   'lon lat time')
      CALL ADD_ATT_REAL(NC_ID, H2OFLX_ID, 'missing_value', msv)

      DIMIDS(1) = lon_Eul_U_DIM
      DIMIDS(2) = lat_Eul_U_DIM
      DIMIDS(3) = depth_UV_DIM
      DIMIDS(4) = time_DIM

      CALL DEF_REAL(NC_ID, 'Eul_U', Eul_U_ID, 4, DIMIDS)
      CALL ADD_ATT_TEXT(NC_ID, Eul_U_ID, 'long_name',
     &   'Eulerian Zonal Velocity')
      CALL ADD_ATT_TEXT(NC_ID, Eul_U_ID, 'units', 'm/s')
      CALL ADD_ATT_TEXT(NC_ID, Eul_U_ID, 'coordinates',
     &   'lon_Eul_U lat_Eul_U depth_UV time')
      CALL ADD_ATT_TEXT(NC_ID, Eul_U_ID, 'vector_angle',
     &   'vector_angle_Eul_U')
      CALL ADD_ATT_REAL(NC_ID, Eul_U_ID, 'missing_value', msv)

      DIMIDS(1) = lon_Eul_U_DIM
      DIMIDS(2) = lat_Eul_U_DIM
      DIMIDS(3) = time_DIM

      CALL DEF_REAL(NC_ID, 'WS_x', WS_x_ID, 3, DIMIDS)
      CALL ADD_ATT_TEXT(NC_ID, WS_x_ID, 'long_name',
     &   'Zonal Wind Stress')
      CALL ADD_ATT_TEXT(NC_ID, WS_x_ID, 'units', 'N/m^2')
      CALL ADD_ATT_TEXT(NC_ID, WS_x_ID, 'coordinates',
     &   'lon_Eul_U lat_Eul_U time')
      CALL ADD_ATT_TEXT(NC_ID, WS_x_ID, 'vector_angle',
     &   'vector_angle_Eul_U')
      CALL ADD_ATT_REAL(NC_ID, WS_x_ID, 'missing_value', msv)

      DIMIDS(1) = lon_Eul_V_DIM
      DIMIDS(2) = lat_Eul_V_DIM
      DIMIDS(3) = depth_UV_DIM
      DIMIDS(4) = time_DIM

      CALL DEF_REAL(NC_ID, 'Eul_V', Eul_V_ID, 4, DIMIDS)
      CALL ADD_ATT_TEXT(NC_ID, Eul_V_ID, 'long_name',
     &   'Eulerian Meridional Velocity')
      CALL ADD_ATT_TEXT(NC_ID, Eul_V_ID, 'units', 'm/s')
      CALL ADD_ATT_TEXT(NC_ID, Eul_V_ID, 'coordinates',
     &   'lon_Eul_V lat_Eul_V depth_UV time')
      CALL ADD_ATT_TEXT(NC_ID, Eul_V_ID, 'vector_angle',
     &   'vector_angle_Eul_V')
      CALL ADD_ATT_REAL(NC_ID, Eul_V_ID, 'missing_value', msv)

      DIMIDS(1) = lon_Eul_V_DIM
      DIMIDS(2) = lat_Eul_V_DIM
      DIMIDS(3) = time_DIM

      CALL DEF_REAL(NC_ID, 'WS_y', WS_y_ID, 3, DIMIDS)
      CALL ADD_ATT_TEXT(NC_ID, WS_y_ID, 'long_name',
     &   'Meridional Wind Stress')
      CALL ADD_ATT_TEXT(NC_ID, WS_y_ID, 'units', 'N/m^2')
      CALL ADD_ATT_TEXT(NC_ID, WS_y_ID, 'coordinates',
     &   'lon_Eul_V lat_Eul_V time')
      CALL ADD_ATT_TEXT(NC_ID, WS_y_ID, 'vector_angle',
     &   'vector_angle_Eul_V')
      CALL ADD_ATT_REAL(NC_ID, WS_y_ID, 'missing_value', msv)

      IF (has_eddy) THEN

         DIMIDS(1) = lon_Eddy_U_DIM
         DIMIDS(2) = lat_Eddy_U_DIM
         DIMIDS(3) = depth_UV_DIM
         DIMIDS(4) = time_DIM

         CALL DEF_REAL(NC_ID, 'Eddy_U', Eddy_U_ID, 4, DIMIDS)
         CALL ADD_ATT_TEXT(NC_ID, Eddy_U_ID, 'long_name',
     &      'Eddy Induced Zonal Velocity')
         CALL ADD_ATT_TEXT(NC_ID, Eddy_U_ID, 'units', 'm/s')
         CALL ADD_ATT_TEXT(NC_ID, Eddy_U_ID, 'coordinates',
     &      'lon_Eddy_U lat_Eddy_U depth_UV time')
         CALL ADD_ATT_TEXT(NC_ID, Eddy_U_ID, 'vector_angle',
     &      'vector_angle_Eddy_U')
         CALL ADD_ATT_REAL(NC_ID, Eddy_U_ID, 'missing_value', msv)

         DIMIDS(1) = lon_Eddy_V_DIM
         DIMIDS(2) = lat_Eddy_V_DIM
         DIMIDS(3) = depth_UV_DIM
         DIMIDS(4) = time_DIM

         CALL DEF_REAL(NC_ID, 'Eddy_V', Eddy_V_ID, 4, DIMIDS)
         CALL ADD_ATT_TEXT(NC_ID, Eddy_V_ID, 'long_name',
     &      'Eddy Induced Meridional Velocity')
         CALL ADD_ATT_TEXT(NC_ID, Eddy_V_ID, 'units', 'm/s')
         CALL ADD_ATT_TEXT(NC_ID, Eddy_V_ID, 'coordinates',
     &      'lon_Eddy_V lat_Eddy_V depth_UV time')
         CALL ADD_ATT_TEXT(NC_ID, Eddy_V_ID, 'vector_angle',
     &      'vector_angle_Eddy_V')
         CALL ADD_ATT_REAL(NC_ID, Eddy_V_ID, 'missing_value', msv)

      ENDIF ! IF (has_eddy)

!***********************************************************************
! Define global attributes
!***********************************************************************
      CALL ADD_ATT_TEXT(NC_ID, NF_GLOBAL, 'output_routine',
     &   '$RCSfile: write_nc_phys.F,v $, $Revision: 1.1.1.1 $')

      CALL ADD_ATT_TEXT(NC_ID, NF_GLOBAL, 'file_name', filename)

      CALL ADD_ATT_TEXT(NC_ID, NF_GLOBAL, 'trx_associate_file',
     &   group//'_grid.nc')
      CALL ADD_ATT_TEXT(NC_ID, NF_GLOBAL, 'vel_associate_file',
     &   group//'_vel_grid.nc')

      CALL ADD_ATT_TEXT(NC_ID, NF_GLOBAL, 'project', 'OCMIP')

      CALL ADD_ATT_TEXT(NC_ID, NF_GLOBAL, 'institution', group)

      CALL ADD_ATT_TEXT(NC_ID, NF_GLOBAL, 'production', production)

      IF (has_eddy) THEN
         CALL ADD_ATT_TEXT(NC_ID, NF_GLOBAL, 'has_eddy', 'TRUE')
      ELSE
         CALL ADD_ATT_TEXT(NC_ID, NF_GLOBAL, 'has_eddy', 'FALSE')
      ENDIF

!***********************************************************************
! End of define mode
!***********************************************************************
      STATUS = NF_ENDDEF(NC_ID)
      IF (STATUS .NE. NF_NOERR) THEN
         PRINT *, 'NetCDF error ending define mode'
         PRINT *, NF_STRERROR(STATUS)
         STOP 'Stopped'
      ENDIF

!***********************************************************************
! Write data to file
!***********************************************************************
      DO I = 1, 12
         IF (I .EQ. 1) THEN
            bounds_TIME(I,1) = 0.0
         ELSE
            bounds_TIME(I,1) = bounds_TIME(I-1,2)
         ENDIF
         TIME(I) = bounds_TIME(I,1) + 0.5 * secs_per_month(i)
         bounds_TIME(I,2) = bounds_TIME(I,1) + secs_per_month(i)
      ENDDO

      CALL PUT_REAL(NC_ID, time_ID, TIME, 'time')
      CALL PUT_REAL(NC_ID, bounds_time_ID, bounds_TIME, 'bounds_TIME')

      CALL PUT_REAL(NC_ID, PT_ID, PT, 'PT')
      CALL PUT_REAL(NC_ID, SAL_ID, SAL, 'SAL')

      CALL PUT_REAL(NC_ID, TFLX_ID, TFLX, 'TFLX')
      CALL PUT_REAL(NC_ID, H2OFLX_ID, H2OFLX, 'H2OFLX')

      CALL PUT_REAL(NC_ID, Eul_U_ID, Eul_U, 'Eul_U')
      CALL PUT_REAL(NC_ID, WS_x_ID, WS_x, 'WS_x')

      CALL PUT_REAL(NC_ID, Eul_V_ID, Eul_V, 'Eul_V')
      CALL PUT_REAL(NC_ID, WS_y_ID, WS_y, 'WS_y')

      IF (has_eddy) THEN
         CALL PUT_REAL(NC_ID, Eddy_U_ID, Eddy_U, 'Eddy_U')
         CALL PUT_REAL(NC_ID, Eddy_V_ID, Eddy_V, 'Eddy_V')
      ENDIF ! IF (has_eddy)

!***********************************************************************
! Close the NetCDF file
!***********************************************************************

      STATUS = NF_CLOSE(NC_ID)
      IF (STATUS .NE. NF_NOERR) THEN
         PRINT *, 'NetCDF error closing ', filename
         PRINT *, NF_STRERROR(STATUS)
         STOP 'Stopped'
      ENDIF

      END