MODULE outncfmd
USE netcdf
USE varsmd
IMPLICIT NONE

  CONTAINS
  SUBROUTINE outncf
    INTEGER :: fid
    INTEGER :: x_dimid, y_dimid, z_dimid
    INTEGER :: z_varid, zz_varid
    INTEGER :: dx_varid, dy_varid
    INTEGER :: east_u_varid, east_v_varid
    INTEGER :: east_e_varid, east_c_varid
    INTEGER :: north_u_varid, north_v_varid
    INTEGER :: north_e_varid, north_c_varid
    INTEGER :: rot_varid, h_varid, fsm_varid
    INTEGER :: dum_varid, dvm_varid
    INTEGER :: t_varid, s_varid
    INTEGER :: u_varid, v_varid, el_varid

    WRITE(gridfile,'(a,''.grid.nc'')')       TRIM(ADJUSTL(runid))
    WRITE(tsinfile,'(a,''.ts_initial.nc'')') TRIM(ADJUSTL(runid)) 
    WRITE(uvinfile,'(a,''.uv_initial.nc'')') TRIM(ADJUSTL(runid))
    WRITE(elinfile,'(a,''.el_initial.nc'')') TRIM(ADJUSTL(runid))
    WRITE(tsclimfile,'(a,''.ts_clim.nc'')')  TRIM(ADJUSTL(runid))

! Create MPIPOM-TC input grid file

    CALL check(NF90_CREATE(gridfile,NF90_CLOBBER,fid))
    CALL check(NF90_PUT_ATT(fid,NF90_GLOBAL,"title",TRIM(ADJUSTL(runid))))
    CALL check(NF90_PUT_ATT(fid,NF90_GLOBAL,"description", &
                                         "POM grid file"))
    CALL check(NF90_DEF_DIM(fid,"x",im,x_dimid))
    CALL check(NF90_DEF_DIM(fid,"y",jm,y_dimid))
    CALL check(NF90_DEF_DIM(fid,"z",kb,z_dimid))
    CALL check(NF90_DEF_VAR(fid,"z",NF90_DOUBLE,z_dimid,z_varid))
    CALL check(NF90_PUT_ATT(fid,z_varid,"long_name", &
                                     "sigma of cell face"))
    CALL check(NF90_PUT_ATT(fid,z_varid,"units","sigma_level"))
    CALL check(NF90_PUT_ATT(fid,z_varid,"standard_name", &
                                     "ocean_sigma_coordinate"))
    CALL check(NF90_PUT_ATT(fid,z_varid,"formula_terms", &
                                     "sigma: z eta: elb depth: h"))
    CALL check(NF90_DEF_VAR(fid,"zz",NF90_DOUBLE,z_dimid,zz_varid))
    CALL check(NF90_PUT_ATT(fid,zz_varid,"long_name", &
                                      "sigma of cell centre"))
    CALL check(NF90_PUT_ATT(fid,zz_varid,"units","sigma_level"))
    CALL check(NF90_PUT_ATT(fid,zz_varid,"standard_name", &
                                      "ocean_sigma_coordinate"))
    CALL check(NF90_PUT_ATT(fid,zz_varid,"formula_terms", &
                                      "sigma: zz eta: elb depth: h"))
    CALL check(NF90_DEF_VAR(fid,"dx",NF90_DOUBLE,(/x_dimid,y_dimid/), &
                                                       dx_varid))
    CALL check(NF90_PUT_ATT(fid,dx_varid,"long_name", &
                                      "grid increment in x"))
    CALL check(NF90_PUT_ATT(fid,dx_varid,"units","metre"))
    CALL check(NF90_PUT_ATT(fid,dx_varid,"coords","east_e north_e"))
    CALL check(NF90_DEF_VAR(fid,"dy",NF90_DOUBLE,(/x_dimid,y_dimid/), &
                                                       dy_varid))
    CALL check(NF90_PUT_ATT(fid,dy_varid,"long_name", &
                                      "grid increment in y"))
    CALL check(NF90_PUT_ATT(fid,dy_varid,"units","metre"))
    CALL check(NF90_PUT_ATT(fid,dy_varid,"coords","east_e north_e"))
    CALL check(NF90_DEF_VAR(fid,"east_u",NF90_DOUBLE, &
                         (/x_dimid,y_dimid/),east_u_varid))
    CALL check(NF90_PUT_ATT(fid,east_u_varid,"long_name", &
                                          "easting of u-points"))
    CALL check(NF90_PUT_ATT(fid,east_u_varid,"units","degree"))
    CALL check(NF90_PUT_ATT(fid,east_u_varid,"coords","east_u north_u"))
    CALL check(NF90_DEF_VAR(fid,"east_v",NF90_DOUBLE, &
                         (/x_dimid,y_dimid/),east_v_varid))
    CALL check(NF90_PUT_ATT(fid,east_v_varid,"long_name", &
                                          "easting of v-points"))
    CALL check(NF90_PUT_ATT(fid,east_v_varid,"units","degree"))
    CALL check(NF90_PUT_ATT(fid,east_v_varid,"coords","east_v north_v"))
    CALL check(NF90_DEF_VAR(fid,"east_e",NF90_DOUBLE, &
                         (/x_dimid,y_dimid/),east_e_varid))
    CALL check(NF90_PUT_ATT(fid,east_e_varid,"long_name", &
                             "easting of elevation points"))
    CALL check(NF90_PUT_ATT(fid,east_e_varid,"units","degree"))
    CALL check(NF90_PUT_ATT(fid,east_e_varid,"coords","east_e north_e"))
    CALL check(NF90_DEF_VAR(fid,"east_c",NF90_DOUBLE, &
                         (/x_dimid,y_dimid/),east_c_varid))
    CALL check(NF90_PUT_ATT(fid,east_c_varid,"long_name", &
                             "easting of cell corners"))
    CALL check(NF90_PUT_ATT(fid,east_c_varid,"units","degree"))
    CALL check(NF90_PUT_ATT(fid,east_c_varid,"coords","east_c north_c"))
    CALL check(NF90_DEF_VAR(fid,"north_u",NF90_DOUBLE, &
                         (/x_dimid,y_dimid/),north_u_varid))
    CALL check(NF90_PUT_ATT(fid,north_u_varid,"long_name", &
                                           "northing of u-points"))
    CALL check(NF90_PUT_ATT(fid,north_u_varid,"units","degree"))
    CALL check(NF90_PUT_ATT(fid,north_u_varid,"coords","east_u north_u"))
    CALL check(NF90_DEF_VAR(fid,"north_v",NF90_DOUBLE, &
                         (/x_dimid,y_dimid/),north_v_varid))
    CALL check(NF90_PUT_ATT(fid,north_v_varid,"long_name", &
                                           "northing of v-points"))
    CALL check(NF90_PUT_ATT(fid,north_v_varid,"units","degree"))
    CALL check(NF90_PUT_ATT(fid,north_v_varid,"coords","east_v north_v"))
    CALL check(NF90_DEF_VAR(fid,"north_e",NF90_DOUBLE, &
                         (/x_dimid,y_dimid/),north_e_varid))
    CALL check(NF90_PUT_ATT(fid,north_e_varid,"long_name", &
                             "northing of elevation points"))
    CALL check(NF90_PUT_ATT(fid,north_e_varid,"units","degree"))
    CALL check(NF90_PUT_ATT(fid,north_e_varid,"coords","east_e north_e"))
    CALL check(NF90_DEF_VAR(fid,"north_c",NF90_DOUBLE, &
                         (/x_dimid,y_dimid/),north_c_varid))
    CALL check(NF90_PUT_ATT(fid,north_c_varid,"long_name", &
                             "northing of cell corners"))
    CALL check(NF90_PUT_ATT(fid,north_c_varid,"units","degree"))
    CALL check(NF90_PUT_ATT(fid,north_c_varid,"coords","east_c north_c"))
    CALL check(NF90_DEF_VAR(fid,"rot",NF90_DOUBLE, &
                         (/x_dimid,y_dimid/),rot_varid))
    CALL check(NF90_PUT_ATT(fid,rot_varid,"long_name", &
                             "rotation angle of x-axis wrt. east"))
    CALL check(NF90_PUT_ATT(fid,rot_varid,"units","degree"))
    CALL check(NF90_PUT_ATT(fid,rot_varid,"coords","east_e north_e"))
    CALL check(NF90_DEF_VAR(fid,"h",NF90_DOUBLE, &
                         (/x_dimid,y_dimid/),h_varid))
    CALL check(NF90_PUT_ATT(fid,h_varid,"long_name", &
                             "undisturbed water depth"))
    CALL check(NF90_PUT_ATT(fid,h_varid,"units","metre"))
    CALL check(NF90_PUT_ATT(fid,h_varid,"coords","east_e north_e"))
    CALL check(NF90_DEF_VAR(fid,"fsm",NF90_DOUBLE, &
                         (/x_dimid,y_dimid/),fsm_varid))
    CALL check(NF90_PUT_ATT(fid,fsm_varid,"long_name", &
                             "free surface mask"))
    CALL check(NF90_PUT_ATT(fid,fsm_varid,"units","dimensionless"))
    CALL check(NF90_PUT_ATT(fid,fsm_varid,"coords","east_e north_e"))
    CALL check(NF90_DEF_VAR(fid,"dum",NF90_DOUBLE, &
                         (/x_dimid,y_dimid/),dum_varid))
    CALL check(NF90_PUT_ATT(fid,dum_varid,"long_name", &
                             "u-velocity mask"))
    CALL check(NF90_PUT_ATT(fid,dum_varid,"units","dimensionless"))
    CALL check(NF90_PUT_ATT(fid,dum_varid,"coords","east_u north_u"))
    CALL check(NF90_DEF_VAR(fid,"dvm",NF90_DOUBLE, &
                         (/x_dimid,y_dimid/),dvm_varid))
    CALL check(NF90_PUT_ATT(fid,dvm_varid,"long_name", &
                             "v-velocity mask"))
    CALL check(NF90_PUT_ATT(fid,dvm_varid,"units","dimensionless"))
    CALL check(NF90_PUT_ATT(fid,dvm_varid,"coords","east_v north_v"))
    CALL check(NF90_ENDDEF(fid))
    CALL check(NF90_PUT_VAR(fid,z_varid,z))
    CALL check(NF90_PUT_VAR(fid,zz_varid,zz))
    CALL check(NF90_PUT_VAR(fid,dx_varid,dx))
    CALL check(NF90_PUT_VAR(fid,dy_varid,dy))
    CALL check(NF90_PUT_VAR(fid,east_u_varid,east_u))
    CALL check(NF90_PUT_VAR(fid,east_v_varid,east_v))
    CALL check(NF90_PUT_VAR(fid,east_c_varid,east_c))
    CALL check(NF90_PUT_VAR(fid,east_e_varid,east_e))
    CALL check(NF90_PUT_VAR(fid,north_u_varid,north_u))
    CALL check(NF90_PUT_VAR(fid,north_v_varid,north_v))
    CALL check(NF90_PUT_VAR(fid,north_c_varid,north_c))
    CALL check(NF90_PUT_VAR(fid,north_e_varid,north_e))
    CALL check(NF90_PUT_VAR(fid,rot_varid,rot))
    CALL check(NF90_PUT_VAR(fid,h_varid,h))
    CALL check(NF90_PUT_VAR(fid,fsm_varid,fsm))
    CALL check(NF90_PUT_VAR(fid,dum_varid,dum))
    CALL check(NF90_PUT_VAR(fid,dvm_varid,dvm))
    CALL check(NF90_CLOSE(fid))

! Create MPIPOM-TC input ts_initial file

    CALL check(NF90_CREATE(tsinfile,NF90_CLOBBER,fid))
    CALL check(NF90_PUT_ATT(fid,NF90_GLOBAL,"title",TRIM(ADJUSTL(runid))))
    CALL check(NF90_PUT_ATT(fid,NF90_GLOBAL,"description", &
                                         "POM initial condition file"))
    CALL check(NF90_DEF_DIM(fid,"x",im,x_dimid))
    CALL check(NF90_DEF_DIM(fid,"y",jm,y_dimid))
    CALL check(NF90_DEF_DIM(fid,"z",nl,z_dimid))
    CALL check(NF90_DEF_VAR(fid,"z",NF90_DOUBLE,z_dimid,z_varid))
    CALL check(NF90_PUT_ATT(fid,z_varid,"long_name","depth level"))
    CALL check(NF90_PUT_ATT(fid,z_varid,"units","m"))
    CALL check(NF90_DEF_VAR(fid,"t",NF90_DOUBLE, &
                         (/x_dimid,y_dimid,z_dimid/),t_varid))
    CALL check(NF90_PUT_ATT(fid,t_varid,"long_name", &
                                     "potential temperature"))
    CALL check(NF90_PUT_ATT(fid,t_varid,"units","K"))
    CALL check(NF90_DEF_VAR(fid,"s",NF90_DOUBLE, &
                         (/x_dimid,y_dimid,z_dimid/),s_varid))
    CALL check(NF90_PUT_ATT(fid,s_varid,"long_name","salinity"))
    CALL check(NF90_PUT_ATT(fid,s_varid,"units","PSS"))
    CALL check(NF90_ENDDEF(fid))
    CALL check(NF90_PUT_VAR(fid,z_varid,zrtof))
    CALL check(NF90_PUT_VAR(fid,t_varid,t))
    CALL check(NF90_PUT_VAR(fid,s_varid,s))
    CALL check(NF90_CLOSE(fid))

! Create MPIPOM-TC input ts_clim file

    CALL check(NF90_CREATE(tsclimfile,NF90_CLOBBER,fid))
    CALL check(NF90_PUT_ATT(fid,NF90_GLOBAL,"title",TRIM(ADJUSTL(runid))))
    CALL check(NF90_PUT_ATT(fid,NF90_GLOBAL,"description", &
                                         "POM initial condition file"))
    CALL check(NF90_DEF_DIM(fid,"x",im,x_dimid))
    CALL check(NF90_DEF_DIM(fid,"y",jm,y_dimid))
    CALL check(NF90_DEF_DIM(fid,"z",nl,z_dimid))
    CALL check(NF90_DEF_VAR(fid,"z",NF90_DOUBLE,z_dimid,z_varid))
    CALL check(NF90_PUT_ATT(fid,z_varid,"long_name","depth level"))
    CALL check(NF90_PUT_ATT(fid,z_varid,"units","m"))
    CALL check(NF90_DEF_VAR(fid,"t",NF90_DOUBLE, &
                         (/x_dimid,y_dimid,z_dimid/),t_varid))
    CALL check(NF90_PUT_ATT(fid,t_varid,"long_name", &
                                     "potential temperature"))
    CALL check(NF90_PUT_ATT(fid,t_varid,"units","K"))
    CALL check(NF90_DEF_VAR(fid,"s",NF90_DOUBLE, &
                         (/x_dimid,y_dimid,z_dimid/),s_varid))
    CALL check(NF90_PUT_ATT(fid,s_varid,"long_name","salinity"))
    CALL check(NF90_PUT_ATT(fid,s_varid,"units","PSS"))
    CALL check(NF90_ENDDEF(fid))
    CALL check(NF90_PUT_VAR(fid,z_varid,zrtof))
    CALL check(NF90_PUT_VAR(fid,t_varid,tclim))
    CALL check(NF90_PUT_VAR(fid,s_varid,sclim))
    CALL check(NF90_CLOSE(fid))

! Create MPIPOM-TC input uv_initial file

    CALL check(NF90_CREATE(uvinfile,NF90_CLOBBER,fid))
    CALL check(NF90_PUT_ATT(fid,NF90_GLOBAL,"title",TRIM(ADJUSTL(runid))))
    CALL check(NF90_PUT_ATT(fid,NF90_GLOBAL,"description", &
                                         "POM initial condition file"))
    CALL check(NF90_DEF_DIM(fid,"x",im,x_dimid))
    CALL check(NF90_DEF_DIM(fid,"y",jm,y_dimid))
    CALL check(NF90_DEF_DIM(fid,"z",nl,z_dimid))
    CALL check(NF90_DEF_VAR(fid,"z",NF90_DOUBLE,z_dimid,z_varid))
    CALL check(NF90_PUT_ATT(fid,z_varid,"long_name","depth level"))
    CALL check(NF90_PUT_ATT(fid,z_varid,"units","m"))
    CALL check(NF90_DEF_VAR(fid,"u",NF90_DOUBLE, &
                         (/x_dimid,y_dimid,z_dimid/),u_varid))
    CALL check(NF90_PUT_ATT(fid,u_varid,"long_name","uvel"))
    CALL check(NF90_PUT_ATT(fid,u_varid,"units","m/s"))
    CALL check(NF90_DEF_VAR(fid,"v",NF90_DOUBLE, &
                         (/x_dimid,y_dimid,z_dimid/),v_varid))
    CALL check(NF90_PUT_ATT(fid,v_varid,"long_name","vvel"))
    CALL check(NF90_PUT_ATT(fid,v_varid,"units","m/s"))
    CALL check(NF90_ENDDEF(fid))
    CALL check(NF90_PUT_VAR(fid,z_varid,zrtof))
    CALL check(NF90_PUT_VAR(fid,u_varid,u))
    CALL check(NF90_PUT_VAR(fid,v_varid,v))
    CALL check(NF90_CLOSE(fid))

! Create MPIPOM-TC input el_initial file

    CALL check(NF90_CREATE(elinfile,NF90_CLOBBER,fid))
    CALL check(NF90_PUT_ATT(fid,NF90_GLOBAL,"title",TRIM(ADJUSTL(runid))))
    CALL check(NF90_PUT_ATT(fid,NF90_GLOBAL,"description", &
                                         "POM initial condition file"))
    CALL check(NF90_DEF_DIM(fid,"x",im,x_dimid))
    CALL check(NF90_DEF_DIM(fid,"y",jm,y_dimid))
    CALL check(NF90_DEF_VAR(fid,"el",NF90_DOUBLE, &
                         (/x_dimid,y_dimid/),el_varid))
    CALL check(NF90_PUT_ATT(fid,el_varid,"long_name","ssh"))
    CALL check(NF90_PUT_ATT(fid,el_varid,"units","m"))
    CALL check(NF90_ENDDEF(fid))
    CALL check(NF90_PUT_VAR(fid,el_varid,el))
    CALL check(NF90_CLOSE(fid))

  END SUBROUTINE outncf

  SUBROUTINE check(index)
    INTEGER, INTENT(IN)                 ::  index
    IF( index /= nf90_noerr) THEN
      PRINT*, nf90_strerror(index)
      STOP
    ENDIF
  END SUBROUTINE check

END MODULE outncfmd