init_3d_ncdf.F90 13.4 KB
Newer Older
kbk's avatar
kbk committed
1
!$Id: init_3d_ncdf.F90,v 1.7 2004-06-15 08:25:57 kbk Exp $
gotm's avatar
gotm committed
2 3 4 5 6 7 8 9 10 11 12 13 14
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: init_3d_ncdf -
!
! !INTERFACE:
   subroutine init_3d_ncdf(fn,title,starttime)
!
! !DESCRIPTION:
!
! !USES:
   use ncdf_3d
15 16 17 18
   use domain, only: ioff,joff,imin,imax,jmin,jmax
#if ! defined(SPHERICAL)
   use domain, only: xc,yc
#endif
gotm's avatar
gotm committed
19 20 21 22 23
   use domain, only: iimin,iimax,jjmin,jjmax,kmax
   use domain, only: grid_type,vert_cord
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
kbk's avatar
kbk committed
24
   character(len=*), intent(in)        :: fn,title,starttime
gotm's avatar
gotm committed
25 26 27 28 29 30 31 32
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
! !REVISION HISTORY:
!
!  $Log: init_3d_ncdf.F90,v $
kbk's avatar
kbk committed
33 34 35 36
!  Revision 1.7  2004-06-15 08:25:57  kbk
!  added supoort for spm - Ruiz
!
!  Revision 1.6  2004/05/04 09:23:51  kbk
37 38 39
!  hydrostatic consistency criteria stored in .3d.nc file
!
!  Revision 1.5  2003/12/16 12:51:04  kbk
40 41 42
!  preparing for proper support for SPM (manuel)
!
!  Revision 1.4  2003/05/09 11:38:26  kbk
43 44 45
!  added proper undef support - based on Adolf Stips patch
!
!  Revision 1.3  2003/04/23 11:53:24  kbk
kbk's avatar
kbk committed
46 47 48
!  save lat/lon info for spherical grid
!
!  Revision 1.2  2003/04/07 12:51:26  kbk
49 50 51 52
!  CURVILINEAR --> defined(SPHERICAL) || defined(CURVILINEAR)
!
!  Revision 1.1.1.1  2002/05/02 14:01:46  gotm
!  recovering after CVS crash
gotm's avatar
gotm committed
53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
!
!  Revision 1.7  2001/10/25 16:16:20  bbh
!  No actual storing of data in init_3d_ncdf.F90 -> save_3d_ncdf.F90
!
!  Revision 1.6  2001/10/23 14:19:20  bbh
!  Stores h if general vertical coordinates
!
!  Revision 1.5  2001/10/23 07:37:17  bbh
!  Saving spm - if calc_spm and save_spm are both true
!
!  Revision 1.4  2001/10/22 08:03:13  bbh
!  Misplaced #else
!
!  Revision 1.3  2001/09/24 14:13:25  bbh
!  xc and yc have changing shape depending on grid_type
!
!  Revision 1.2  2001/09/19 11:20:32  bbh
!  Explicit de-allocates memory when -DFORTRAN90
!
!  Revision 1.1  2001/09/13 14:50:02  bbh
!  Cleaner and smaller NetCDF implementation + better axis support
!
! !LOCAL VARIABLES:
kbk's avatar
kbk committed
76 77 78 79 80 81 82
   integer                   :: i,j,k,err
   character(32)             :: xname,yname,zname,xunits,yunits,zunits
   character(32)             :: gridname,vertname
   integer                   :: scalar(1),axisdim(1),f3_dims(3),f4_dims(4)
   REALTYPE                  :: fv,mv,vr(2)
   REALTYPE                  :: x
   character(len=80)         :: history,ts
gotm's avatar
gotm committed
83 84 85 86 87 88 89 90 91 92 93 94 95 96 97
!EOP
!-------------------------------------------------------------------------
!BOC
   include "netcdf.inc"

   xlen = iimax-iimin+1
   ylen = jjmax-jjmin+1
   zlen = kmax+1

   select case (grid_type)
      case (1)
         xname = 'xax'
         yname = 'yax'
         xunits = 'meters'
         yunits = 'meters'
kbk's avatar
kbk committed
98
         gridname='cartesian'
gotm's avatar
gotm committed
99 100 101 102 103
      case (2)
         xname = 'lon'
         yname = 'lat'
         xunits = 'degrees_east'
         yunits = 'degrees_north'
kbk's avatar
kbk committed
104
         gridname='spherical'
gotm's avatar
gotm committed
105 106 107 108 109
      case (3)
         xname = 'xdim'
         yname = 'ydim'
         xunits = 'curvi-x'
         yunits = 'curvi-y'
kbk's avatar
kbk committed
110
         gridname='curvi-linear'
gotm's avatar
gotm committed
111 112 113 114 115 116 117
      case default
   end select

   select case (vert_cord)
      case (1)
         zname = 'zax'
         zunits = 'sigma_level'
kbk's avatar
kbk committed
118
         vertname='sigma coordinates'
gotm's avatar
gotm committed
119 120 121
      case (2)
         zname = 'zax'
         zunits = 'meters'
kbk's avatar
kbk committed
122
         vertname='z-levels'
gotm's avatar
gotm committed
123 124 125
      case (3)
         zname = 'zax'
         zunits = 'level'
kbk's avatar
kbk committed
126
         vertname='general vertical coordinates'
gotm's avatar
gotm committed
127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145
      case default
   end select

   err = nf_create(fn, NF_CLOBBER, ncid)
   if (err .NE. NF_NOERR) go to 10

!  dimensions
   err = nf_def_dim(ncid,xname,xlen,x_dim)
   if (err .NE. NF_NOERR) go to 10

   err = nf_def_dim(ncid,yname,ylen,y_dim)
   if (err .NE. NF_NOERR) go to 10

   err = nf_def_dim(ncid,zname,zlen,z_dim)
   if (err .NE. NF_NOERR) go to 10

   err = nf_def_dim(ncid,'time',NF_UNLIMITED,time_dim)
   if (err .NE. NF_NOERR) go to 10

kbk's avatar
kbk committed
146 147 148
   f3_dims(3)= time_dim
   f3_dims(2)= y_dim
   f3_dims(1)= x_dim
gotm's avatar
gotm committed
149

kbk's avatar
kbk committed
150 151 152 153
   f4_dims(4)= time_dim
   f4_dims(3)= z_dim
   f4_dims(2)= y_dim
   f4_dims(1)= x_dim
gotm's avatar
gotm committed
154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182

   history = 'Generated by getm, ver. '//RELEASE
   ts = 'seconds since '//starttime

!  info on offset, grid type and vertical coordinates
   scalar(1) = 0
   err = nf_def_var(ncid,'grid_type',NF_INT,0,scalar,grid_type_id)
   if (err .NE. NF_NOERR) go to 10
   call set_attributes(ncid,grid_type_id,units=gridname)

   err = nf_def_var(ncid,'vert_cord',NF_INT,0,scalar,vert_cord_id)
   if (err .NE. NF_NOERR) go to 10
   call set_attributes(ncid,vert_cord_id,units=vertname)

   err = nf_def_var(ncid,'ioff',NF_INT,0,scalar,ioff_id)
   if (err .NE. NF_NOERR) go to 10
   call set_attributes(ncid,ioff_id,long_name='index offset (i)')

   err = nf_def_var(ncid,'joff',NF_INT,0,scalar,joff_id)
   if (err .NE. NF_NOERR) go to 10
   call set_attributes(ncid,joff_id,long_name='index offset (j)')

!  time
   axisdim(1) = time_dim
   err = nf_def_var(ncid,'time',NF_REAL,1,axisdim,time_id)
   if (err .NE. NF_NOERR) go to 10
   call set_attributes(ncid,time_id,units=trim(ts),long_name='time')

!  coordinate variables
kbk's avatar
kbk committed
183 184
   select case (grid_type)
      case (1)
185
#if ! ( defined(SPHERICAL) || defined(CURVILINEAR) )
kbk's avatar
kbk committed
186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207
         axisdim(1) = x_dim
         err = nf_def_var(ncid,xname,NF_REAL,1,axisdim,xc_id)
         if (err .NE. NF_NOERR) go to 10
         call set_attributes(ncid,xc_id,units=xunits)
         axisdim(1) = y_dim
         err = nf_def_var(ncid,yname,NF_REAL,1,axisdim,yc_id)
         if (err .NE. NF_NOERR) go to 10
         call set_attributes(ncid,yc_id,units=yunits)
#endif
      case (2)
#if defined(SPHERICAL)
         axisdim(1) = x_dim
         err = nf_def_var(ncid,xname,NF_REAL,1,axisdim,lonc_id)
         if (err .NE. NF_NOERR) go to 10
         call set_attributes(ncid,lonc_id,units=xunits)
         axisdim(1) = y_dim
         err = nf_def_var(ncid,yname,NF_REAL,1,axisdim,latc_id)
         if (err .NE. NF_NOERR) go to 10
         call set_attributes(ncid,latc_id,units=yunits)
#endif
      case (3)
#if defined(CURVILINEAR)
gotm's avatar
gotm committed
208 209
!  do something about curvi-linear coordinates - define xu,xv,....
#endif
kbk's avatar
kbk committed
210 211 212
      case default
   end select

gotm's avatar
gotm committed
213 214 215 216 217 218 219 220
   axisdim(1) = z_dim
   err = nf_def_var(ncid,zname,NF_REAL,1,axisdim,z_id)
   if (err .NE. NF_NOERR) go to 10
   call set_attributes(ncid,z_id,units=zunits)

!  bathymetry
   err = nf_def_var(ncid,'bathymetry',NF_REAL,2,f3_dims,bathymetry_id)
   if (err .NE. NF_NOERR) go to 10
221 222
   fv = h_missing
   mv = h_missing
gotm's avatar
gotm committed
223 224
   vr(1) = -5.
   vr(2) = 4000.
kbk's avatar
kbk committed
225 226 227
   call set_attributes(ncid,bathymetry_id,  &
                       long_name='bathymetry',units='meters',          &
                       FillValue=fv,missing_value=mv,valid_range=vr)
gotm's avatar
gotm committed
228

229 230 231 232 233 234 235 236 237 238 239
!  hydrostatic consistency criteria
   err = nf_def_var(ncid,'hcc',NF_REAL,3,f4_dims,hcc_id)
   if (err .NE. NF_NOERR) go to 10
   fv = -_ONE_ 
   mv = -_ONE_
   vr(1) = 0.
   vr(2) = 1.
   call set_attributes(ncid,hcc_id,  &
                       long_name='hcc',units=' ',          &
                       FillValue=fv,missing_value=mv,valid_range=vr)

gotm's avatar
gotm committed
240 241 242 243
!  now to the variables
!  elevation
   err = nf_def_var(ncid,'elev',NF_REAL,3,f3_dims,elev_id)
   if (err .NE. NF_NOERR) go to 10
244 245
   fv = elev_missing
   mv = elev_missing
gotm's avatar
gotm committed
246 247 248
   vr(1) = -15.
   vr(2) =  15.
   call set_attributes(ncid,elev_id,long_name='elevation',units='meters', &
kbk's avatar
kbk committed
249
                       FillValue=fv,missing_value=mv,valid_range=vr)
gotm's avatar
gotm committed
250 251 252 253

!  depth integrated zonal velocity
   err = nf_def_var(ncid,'u',NF_REAL,3,f3_dims,u_id)
   if (err .NE. NF_NOERR) go to 10
254 255
   fv = vel_missing
   mv = vel_missing
gotm's avatar
gotm committed
256 257 258
   vr(1) = -3.
   vr(2) =  3.
   call set_attributes(ncid,u_id,long_name='int. zonal vel.',units='m/s', &
kbk's avatar
kbk committed
259
                       FillValue=fv,missing_value=mv,valid_range=vr)
gotm's avatar
gotm committed
260

261
!  depth integrated meridional velocity
gotm's avatar
gotm committed
262 263
   err = nf_def_var(ncid,'v',NF_REAL,3,f3_dims,v_id)
   if (err .NE. NF_NOERR) go to 10
264 265
   fv = vel_missing
   mv = vel_missing
gotm's avatar
gotm committed
266 267 268
   vr(1) = -3.
   vr(2) =  3.
   call set_attributes(ncid,v_id,long_name='int. meridional vel.',units='m/s', &
kbk's avatar
kbk committed
269
                       FillValue=fv,missing_value=mv,valid_range=vr)
gotm's avatar
gotm committed
270 271 272 273 274

   select case (vert_cord)
      case (1)
      case (2)
         STDERR 'store z-levels'
kbk's avatar
kbk committed
275
         stop 'init_3d_ncdf'
gotm's avatar
gotm committed
276
      case (3)
277 278
         fv = h_missing
         mv = h_missing
gotm's avatar
gotm committed
279 280
         err = nf_def_var(ncid,'h',NF_REAL,4,f4_dims,h_id)
         if (err .NE. NF_NOERR) go to 10
kbk's avatar
kbk committed
281 282
         call set_attributes(ncid,h_id,long_name='layer thickness',  &
                             units='meters',FillValue=fv,missing_value=mv)
gotm's avatar
gotm committed
283 284 285 286
      case default
   end select

   if (save_vel) then
287 288
      fv = vel_missing
      mv = vel_missing
gotm's avatar
gotm committed
289 290 291 292 293 294 295
      vr(1) = -3.
      vr(2) =  3.

!     zonal velocity
      err = nf_def_var(ncid,'uu',NF_REAL,4,f4_dims,uu_id)
      if (err .NE. NF_NOERR) go to 10
      call set_attributes(ncid,uu_id,long_name='zonal vel.',units='m/s', &
kbk's avatar
kbk committed
296
                          FillValue=fv,missing_value=mv,valid_range=vr)
gotm's avatar
gotm committed
297 298 299 300 301

!     meridional velocity
      err = nf_def_var(ncid,'vv',NF_REAL,4,f4_dims,vv_id)
      if (err .NE. NF_NOERR) go to 10
      call set_attributes(ncid,vv_id,long_name='meridional vel.',units='m/s', &
kbk's avatar
kbk committed
302
                          FillValue=fv,missing_value=mv,valid_range=vr)
gotm's avatar
gotm committed
303 304 305 306 307

!     vertical velocity
      err = nf_def_var(ncid,'w',NF_REAL,4,f4_dims,w_id)
      if (err .NE. NF_NOERR) go to 10
      call set_attributes(ncid,w_id,long_name='vertical vel.',units='m/s', &
kbk's avatar
kbk committed
308
                          FillValue=fv,missing_value=mv,valid_range=vr)
gotm's avatar
gotm committed
309 310 311 312 313
   end if

   if (save_strho) then

      if (save_s) then
314 315
         fv = salt_missing
         mv = salt_missing
gotm's avatar
gotm committed
316 317 318 319 320
         vr(1) =  0.
         vr(2) = 40.
         err = nf_def_var(ncid,'salt',NF_REAL,4,f4_dims,salt_id)
         if (err .NE. NF_NOERR) go to 10
         call set_attributes(ncid,salt_id,long_name='salinity',units='PSU', &
kbk's avatar
kbk committed
321
                             FillValue=fv,missing_value=mv,valid_range=vr)
gotm's avatar
gotm committed
322 323 324
      end if

      if (save_t) then
325 326
         fv = temp_missing
         mv = temp_missing
gotm's avatar
gotm committed
327 328 329 330 331
         vr(1) =  0.
         vr(2) = 40.
         err = nf_def_var(ncid,'temp',NF_REAL,4,f4_dims,temp_id)
         if (err .NE. NF_NOERR) go to 10
         call set_attributes(ncid,temp_id,long_name='temperature',units='degC',&
kbk's avatar
kbk committed
332
                             FillValue=fv,missing_value=mv,valid_range=vr)
gotm's avatar
gotm committed
333 334 335
      end if

      if (save_rho) then
336 337
         fv = rho_missing
         mv = rho_missing
gotm's avatar
gotm committed
338 339 340 341 342
         vr(1) =  0.
         vr(2) = 30.
         err = nf_def_var(ncid,'sigma_t',NF_REAL,4,f4_dims,sigma_t_id)
         if (err .NE. NF_NOERR) go to 10
         call set_attributes(ncid,sigma_t_id,long_name='sigma_t',units='kg/m3',&
kbk's avatar
kbk committed
343
                             FillValue=fv,missing_value=mv,valid_range=vr)
gotm's avatar
gotm committed
344 345 346 347 348 349
      end if
   end if

   if (save_turb) then

      if (save_tke) then
350 351
         fv = tke_missing
         mv = tke_missing
gotm's avatar
gotm committed
352 353 354 355 356
         vr(1) = 0.
         vr(2) = 0.2
         err = nf_def_var(ncid,'tke',NF_REAL,4,f4_dims,tke_id)
         if (err .NE. NF_NOERR) go to 10
         call set_attributes(ncid,tke_id,long_name='TKE',units='m2/s2', &
kbk's avatar
kbk committed
357
                             FillValue=fv,missing_value=mv,valid_range=vr)
gotm's avatar
gotm committed
358 359 360
      end if

      if (save_num) then
361 362
         fv = num_missing
         mv = num_missing
gotm's avatar
gotm committed
363 364 365 366 367
         vr(1) = 0.
         vr(2) = 0.2
         err = nf_def_var(ncid,'num',NF_REAL,4,f4_dims,num_id)
         if (err .NE. NF_NOERR) go to 10
         call set_attributes(ncid,num_id,long_name='viscosity',units='m2/s', &
kbk's avatar
kbk committed
368
                             FillValue=fv,missing_value=mv,valid_range=vr)
gotm's avatar
gotm committed
369 370 371
      end if

!      if (save_nuh) then
372 373
         fv = nuh_missing
         mv = nuh_missing
gotm's avatar
gotm committed
374 375 376 377 378
         vr(1) = 0.
         vr(2) = 0.2
         err = nf_def_var(ncid,'nuh',NF_REAL,4,f4_dims,nuh_id)
         if (err .NE. NF_NOERR) go to 10
         call set_attributes(ncid,nuh_id,long_name='diffusivity',units='m2/s', &
kbk's avatar
kbk committed
379
                             FillValue=fv,missing_value=mv,valid_range=vr)
gotm's avatar
gotm committed
380 381 382
!      end if

      if (save_eps) then
383 384
         fv = eps_missing
         mv = eps_missing
gotm's avatar
gotm committed
385 386 387 388 389
         vr(1) = 0.
         vr(2) = 0.2
         err = nf_def_var(ncid,'diss',NF_REAL,4,f4_dims,eps_id)
         if (err .NE. NF_NOERR) go to 10
         call set_attributes(ncid,eps_id,long_name='dissipation',units='m2/s3',&
kbk's avatar
kbk committed
390
                             FillValue=fv,missing_value=mv,valid_range=vr)
gotm's avatar
gotm committed
391 392 393
      end if
   end if

394
#ifdef SPM
gotm's avatar
gotm committed
395
   if (save_spm) then
396 397
      fv = spm_missing
      mv = spm_missing
kbk's avatar
kbk committed
398 399 400 401 402 403 404
      err = nf_def_var(ncid,'spm_pool',NF_REAL,3,f3_dims,spmpool_id) 
      if (err .NE. NF_NOERR) go to 10
      vr(1) = 0.
      vr(2) = 10.
      call set_attributes(ncid,spmpool_id,long_name='bottom spm pool', &
                          units='kg/m2', & 
                          FillValue=fv,missing_value=mv,valid_range=vr)
gotm's avatar
gotm committed
405 406 407 408
      vr(1) =  0.
      vr(2) = 30.
      err = nf_def_var(ncid,'spm',NF_REAL,4,f4_dims,spm_id)
      if (err .NE. NF_NOERR) go to 10
kbk's avatar
kbk committed
409
      call set_attributes(ncid,spm_id,  &
kbk's avatar
kbk committed
410 411 412
                          long_name='suspended particulate matter', &
                          units='kg/m3', &
                          FillValue=fv,missing_value=mv,valid_range=vr)
gotm's avatar
gotm committed
413
   end if
414
#endif
kbk's avatar
kbk committed
415

gotm's avatar
gotm committed
416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436
!  globals
   err = nf_put_att_text(ncid,NF_GLOBAL,'title',LEN_TRIM(title),title)
   if (err .NE. NF_NOERR) go to 10

   err = nf_put_att_text(ncid,NF_GLOBAL,'history',LEN_TRIM(history),history)
   if (err .NE. NF_NOERR) go to 10

   ! leave define mode
   err = nf_enddef(ncid)
   if (err .NE. NF_NOERR) go to 10

   return

   10 FATAL 'init_3d_ncdf: ',nf_strerror(err)
   stop 'init_3d_ncdf'
   end subroutine init_3d_ncdf
!EOC

!-----------------------------------------------------------------------
! Copyright (C) 2001 - Hans Burchard and Karsten Bolding (BBH)         !
!-----------------------------------------------------------------------