← Back to branch summary

~fluidity-core/fluidity/sea-ice-branch

« back to all changes in this revision

Viewing changes to parameterisation/Equation_of_State.F90

  • Committer: Simon Mouradian
  • Date: 2012-10-19 10:35:59 UTC
  • mfrom: (3520.32.371 fluidity)
  • Revision ID: simon.mouradian06@imperial.ac.uk-20121019103559-y36qa47phc69q8sc
merge from trunk

Show diffs side-by-side

added added

removed removed

Lines of Context:
parameterisation/Equation_of_State.F90
32
32
  use fields
33
33
  use state_module
34
34
  use global_parameters, only: OPTION_PATH_LEN
 
35
  use sediment, only: get_n_sediment_fields, get_sediment_item
35
36
  use spud
36
 
  use sediment, only: get_sediment_name, get_nSediments
37
37
  implicit none
38
38
  
39
39
  private
54
54
    type(scalar_field), pointer:: T, S, oldT, oldS, topdis
55
55
    type(scalar_field) DeltaT, DeltaS, remapT, remapS, fluidconcentration,&
56
56
         & sedimentdensity
57
 
    character(len=OPTION_PATH_LEN) option_path, dep_option_path, class_name, sfield_name
 
57
    character(len=OPTION_PATH_LEN) option_path, dep_option_path, sediment_field_name, class_name, sfield_name
58
58
    logical, dimension(:), allocatable:: done
59
59
    logical include_depth_below
60
60
    real T0, S0, gamma, rho_0, salt, temp, dist, dens, theta
61
61
    integer, dimension(:), pointer:: density_nodes
62
 
    integer ele, i, node, nSediments, f
 
62
    integer ele, i, node, n_sediment_fields, f
63
63
    
64
64
    ewrite(1,*) 'In calculate_perturbation_density'
65
65
    
205
205
       call allocate(deltaS, density%mesh, "DeltaS")
206
206
       call allocate(remapS, density%mesh, "RemapS")
207
207
       call allocate(sedimentdensity, density%mesh, "SedimentDensity")
208
 
       ! fluidconcentration is 1-sedimentconcentration.
209
 
       call allocate(fluidconcentration, density%mesh, &
210
 
            "FluidConcentration")
211
208
       call zero(sedimentdensity)
212
 
       call set(fluidconcentration, 1.0)
213
 
       nSediments=get_nSediments()
214
 
 
215
 
       do i=1,nSediments
216
 
 
217
 
          class_name = get_sediment_name(i)
218
 
 
219
 
          S=>extract_scalar_field(state,trim(class_name))
220
 
          option_path = S%option_path
221
 
 
222
 
          call get_option(trim(option_path)//'/density', gamma)
223
 
          ! Note that 1000.0 is a hack. We actually need to properly
224
 
          ! account for the reference density of seawater.
225
 
          gamma=rho_0*(gamma/1000.0) - rho_0
 
209
 
 
210
       n_sediment_fields = get_n_sediment_fields()
 
211
 
 
212
       do i=1,n_sediment_fields
 
213
       
 
214
          call get_sediment_item(state, i, S)
 
215
 
 
216
          call get_sediment_item(state, i, 'submerged_specific_gravity', gamma)
 
217
          
 
218
          gamma = gamma * rho_0
226
219
 
227
220
          oldS => extract_scalar_field(state, &
228
 
               "Old"//trim(class_name))
 
221
               "Old"//trim(S%name))
229
222
          
230
223
          ! deltaS=theta*S+(1-theta)*oldS-S0
231
224
          call remap_field(S, remapS)
236
229
          call addto(deltaS, remapS, 1.0-theta)
237
230
          ! density=density+gamma*deltaS
238
231
          call addto(sedimentdensity, deltaS, scale=gamma)
239
 
          call addto(fluidconcentration, deltaS, scale=-1.0)
240
232
 
241
233
       end do
242
234
       
243
 
       call scale(density,fluidconcentration)
244
235
       call addto(density,sedimentdensity)
245
236
 
246
237
       call deallocate(deltaS)
247
238
       call deallocate(remapS)
248
 
       call deallocate(fluidconcentration)
249
239
       call deallocate(sedimentdensity)
250
240
       
251
241
    end if