~ubuntu-branches/ubuntu/vivid/cdftools/vivid

« back to all changes in this revision

Viewing changes to cdfspice.f90

  • Committer: Package Import Robot
  • Author(s): Alastair McKinstry
  • Date: 2013-11-14 19:04:43 UTC
  • Revision ID: package-import@ubuntu.com-20131114190443-ymhovvnzvr5kd02l
Tags: upstream-3.0
ImportĀ upstreamĀ versionĀ 3.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
PROGRAM cdfspice
 
2
  !!======================================================================
 
3
  !!                     ***  PROGRAM  cdfspice  ***
 
4
  !!=====================================================================
 
5
  !!  ** Purpose : Compute spiciness 3D field from gridT file
 
6
  !!               Store the results on a 'similar' cdf file.
 
7
  !!
 
8
  !!  ** Method  :  spiciness = sum(i=0,5)[sum(j=0,4)[b(i,j)*theta^i*(s-35)^j]]
 
9
  !!                   with:  b     -> coefficients
 
10
  !!                          theta -> potential temperature
 
11
  !!                          s     -> salinity
 
12
  !!
 
13
  !!  **  Example:
 
14
  !!       spice(15,33)=   0.5445863      0.544586321373410  calcul en double
 
15
  !!       spice(15,33)=   0.5445864      (calcul en simple precision)
 
16
  !!
 
17
  !!  ** References : Flament (2002) "A state variable for characterizing 
 
18
  !!              water masses and their diffusive stability: spiciness."
 
19
  !!              Progress in Oceanography Volume 54, 2002, Pages 493-501.
 
20
  !!
 
21
  !! History : 2.1  : 03/2010  : C.O. Dufour  : Original code
 
22
  !!           3.0  : 01/2011  : J.M. Molines : Doctor norm + Lic.
 
23
  !!----------------------------------------------------------------------
 
24
  USE cdfio
 
25
  USE modcdfnames
 
26
  !!----------------------------------------------------------------------
 
27
  !! CDFTOOLS_3.0 , MEOM 2011
 
28
  !! $Id$
 
29
  !! Copyright (c) 2011, J.-M. Molines
 
30
  !! Software governed by the CeCILL licence (Licence/CDFTOOLSCeCILL.txt)
 
31
  !!----------------------------------------------------------------------
 
32
  IMPLICIT NONE
 
33
 
 
34
  INTEGER(KIND=4)                           :: ji, jj, jk, jt     ! dummy loop index
 
35
  INTEGER(KIND=4)                           :: ierr               ! error status
 
36
  INTEGER(KIND=4)                           :: narg, iargc        ! browse command line
 
37
  INTEGER(KIND=4)                           :: npiglo, npjglo     ! size of the domain
 
38
  INTEGER(KIND=4)                           :: npk, npt           ! size of the domain
 
39
  INTEGER(KIND=4)                           :: ncout              ! ncid of output file
 
40
  INTEGER(KIND=4), DIMENSION(1)             :: ipk, id_varout     ! level and  varid's
 
41
 
 
42
  REAL(KIND=4)                              :: zspval             ! missing value
 
43
  REAL(KIND=4), DIMENSION(:),   ALLOCATABLE :: tim                ! time counter
 
44
 
 
45
  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: dtemp              ! temperature
 
46
  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: dtempt             ! temperature
 
47
  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: dsal               ! salinity
 
48
  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: dsalt              ! salinity
 
49
  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: dsalref            ! reference salinity
 
50
  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: dspi               ! spiceness
 
51
  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: dmask              ! 2D mask at current level
 
52
  REAL(KIND=8), DIMENSION(6,5)              :: dbet               ! coefficients of spiciness formula
 
53
 
 
54
  CHARACTER(LEN=256)                        :: cf_tfil            ! input filename
 
55
  CHARACTER(LEN=256)                        :: cf_out='spice.nc'  ! output file name
 
56
 
 
57
  TYPE (variable), DIMENSION(1)             :: stypvar            ! structure for attributes
 
58
  !!----------------------------------------------------------------------
 
59
  CALL ReadCdfNames()
 
60
 
 
61
  narg = iargc()
 
62
  IF ( narg == 0 ) THEN
 
63
     PRINT *,' usage : cdfspice T-file '
 
64
     PRINT *,'      '
 
65
     PRINT *,'     PURPOSE :'
 
66
     PRINT *,'       Compute the spiceness corresponding to temperatures and salinities'
 
67
     PRINT *,'       given in the input file.' 
 
68
     PRINT *,'      '
 
69
     PRINT *,'       spiciness = sum(i=0,5)[sum(j=0,4)[b(i,j)*theta^i*(s-35)^j]]'
 
70
     PRINT *,'                 with:  b     -> coefficients'
 
71
     PRINT *,'                        theta -> potential temperature'
 
72
     PRINT *,'                        s     -> salinity'
 
73
     PRINT *,'      '
 
74
     PRINT *,'     ARGUMENTS :'
 
75
     PRINT *,'       T-file : netcdf file with temperature and salinity (gridT)' 
 
76
     PRINT *,'      '
 
77
     PRINT *,'     REQUIRED FILES :'
 
78
     PRINT *,'       none'
 
79
     PRINT *,'      '
 
80
     PRINT *,'     OUTPUT : '
 
81
     PRINT *,'       netcdf file : ', TRIM(cf_out) 
 
82
     PRINT *,'         variables : vospice'
 
83
     PRINT *,'      '
 
84
     PRINT *,'     REFERENCE :'
 
85
     PRINT *,'       Flament (2002) "A state variable for characterizing '
 
86
     PRINT *,'             water masses and their diffusive stability: spiciness."'
 
87
     PRINT *,'             Progress in Oceanography Volume 54, 2002, Pages 493-501.'
 
88
     STOP
 
89
  ENDIF
 
90
  IF ( narg == 0 ) THEN
 
91
     PRINT *,'usage : cdfspice  gridT '
 
92
     PRINT *,'    Output on spice.nc, variable vospice'
 
93
     STOP
 
94
  ENDIF
 
95
 
 
96
  CALL getarg (1, cf_tfil)
 
97
 
 
98
  IF ( chkfile(cf_tfil) ) STOP ! missing files
 
99
 
 
100
  npiglo = getdim (cf_tfil,cn_x)
 
101
  npjglo = getdim (cf_tfil,cn_y)
 
102
  npk    = getdim (cf_tfil,cn_z)
 
103
  npt    = getdim (cf_tfil,cn_t)
 
104
 
 
105
  ipk(:)                       = npk 
 
106
  stypvar(1)%cname             = 'vospice'
 
107
  stypvar(1)%cunits            = 'kg/m3'
 
108
  stypvar(1)%rmissing_value    = 0.
 
109
  stypvar(1)%valid_min         = -300.
 
110
  stypvar(1)%valid_max         = 300.
 
111
  stypvar(1)%clong_name        = 'spiciness'
 
112
  stypvar(1)%cshort_name       = 'vospice'
 
113
  stypvar(1)%conline_operation = 'N/A'
 
114
  stypvar(1)%caxis             = 'TZYX'
 
115
 
 
116
  PRINT *, 'npiglo = ', npiglo
 
117
  PRINT *, 'npjglo = ', npjglo
 
118
  PRINT *, 'npk    = ', npk
 
119
  PRINT *, 'npt    = ', npt
 
120
 
 
121
  ALLOCATE (dtemp(npiglo,npjglo), dsal (npiglo,npjglo) )
 
122
  ALLOCATE (dspi( npiglo,npjglo), dmask(npiglo,npjglo) )
 
123
  ALLOCATE (dtempt(npiglo,npjglo), dsalt(npiglo,npjglo))
 
124
  ALLOCATE (dsalref(npiglo,npjglo))
 
125
  ALLOCATE (tim(npt))
 
126
 
 
127
  ! create output fileset
 
128
  ncout = create      (cf_out, cf_tfil, npiglo, npjglo, npk       )
 
129
  ierr  = createvar   (ncout,  stypvar, 1,      ipk,    id_varout )
 
130
  ierr  = putheadervar(ncout,  cf_tfil, npiglo, npjglo, npk       )
 
131
 
 
132
  tim  = getvar1d(cf_tfil, cn_vtimec, npt     )
 
133
  ierr = putvar1d(ncout,   tim,       npt, 'T')
 
134
 
 
135
  zspval = getatt(cf_tfil, cn_vosaline, 'missing_value')
 
136
 
 
137
  ! Define coefficients to compute spiciness (R*8)
 
138
  dbet(1,1) = 0           ; dbet(1,2) = 7.7442d-01  ; dbet(1,3) = -5.85d-03   ; dbet(1,4) = -9.84d-04   ; dbet(1,5) = -2.06d-04
 
139
  dbet(2,1) = 5.1655d-02  ; dbet(2,2) = 2.034d-03   ; dbet(2,3) = -2.742d-04  ; dbet(2,4) = -8.5d-06    ; dbet(2,5) = 1.36d-05
 
140
  dbet(3,1) = 6.64783d-03 ; dbet(3,2) = -2.4681d-04 ; dbet(3,3) = -1.428d-05  ; dbet(3,4) = 3.337d-05   ; dbet(3,5) = 7.894d-06
 
141
  dbet(4,1) = -5.4023d-05 ; dbet(4,2) = 7.326d-06   ; dbet(4,3) = 7.0036d-06  ; dbet(4,4) = -3.0412d-06 ; dbet(4,5) = -1.0853d-06
 
142
  dbet(5,1) = 3.949d-07   ; dbet(5,2) = -3.029d-08  ; dbet(5,3) = -3.8209d-07 ; dbet(5,4) = 1.0012d-07  ; dbet(5,5) = 4.7133d-08
 
143
  dbet(6,1) = -6.36d-10   ; dbet(6,2) = -1.309d-09  ; dbet(6,3) = 6.048d-09   ; dbet(6,4) = -1.1409d-09 ; dbet(6,5) = -6.676d-10
 
144
 
 
145
  ! Compute spiciness
 
146
  DO jt=1,npt
 
147
     PRINT *,' TIME = ', jt, tim(jt)/86400.,' days'
 
148
     DO jk = 1, npk
 
149
        dmask(:,:) = 1.
 
150
 
 
151
        dtemp(:,:) = getvar(cf_tfil, cn_votemper, jk, npiglo, npjglo, ktime=jt)
 
152
        dsal( :,:) = getvar(cf_tfil, cn_vosaline, jk, npiglo, npjglo, ktime=jt)
 
153
 
 
154
        WHERE(dsal == zspval ) dmask = 0
 
155
 
 
156
        ! spiciness at time jt, at level jk  
 
157
        dspi(:,:)    = 0.d0
 
158
        dsalref(:,:) = dsal(:,:) - 35.d0
 
159
        dtempt(:,:)  = 1.d0
 
160
        DO ji=1,6
 
161
           dsalt(:,:) = 1.d0
 
162
           DO jj=1,5
 
163
              dspi( :,:) = dspi (:,:) + dbet   (ji,jj) * dtempt(:,:) * dsalt(:,:)
 
164
              dsalt(:,:) = dsalt(:,:) * dsalref( :,: )
 
165
           END DO
 
166
           dtempt(:,:) = dtempt(:,:) * dtemp(:,:)     
 
167
        END DO
 
168
 
 
169
        ierr = putvar(ncout, id_varout(1), REAL(dspi*dmask), jk, npiglo, npjglo, ktime=jt)
 
170
 
 
171
     END DO  ! loop to next level
 
172
  END DO  ! next time frame
 
173
 
 
174
  ierr = closeout(ncout)
 
175
 
 
176
END PROGRAM cdfspice