← Back to branch summary

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

« back to all changes in this revision

Viewing changes to main/Fluids.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:
main/Fluids.F90
109
109
  use sea_ice_module
110
110
#endif
111
111
  use detector_parallel, only: sync_detector_coordinates, deallocate_detector_list_array
112
 
 
 
112
  use momentum_diagnostic_fields, only: calculate_densities
113
113
 
114
114
  implicit none
115
115
 
168
168
    INTEGER :: ss,ph
169
169
    LOGICAL :: have_solids
170
170
 
171
 
    !     Turbulence modelling - JBull 24-05-11
172
 
    LOGICAL :: have_k_epsilon
173
 
    character(len=OPTION_PATH_LEN) :: keps_option_path
 
171
    ! An array of submaterials of the current phase in state(istate).
 
172
    ! Needed for k-epsilon VelocityBuoyancyDensity calculation line:~630
 
173
    ! S Parkinson 31-08-12
 
174
    type(state_type), dimension(:), pointer :: submaterials     
174
175
 
175
176
    ! Pointers for scalars and velocity fields
176
177
    type(scalar_field), pointer :: sfield
201
202
    call initialise_qmesh
202
203
    call initialise_write_state
203
204
 
204
 
 
205
 
    ! Initialise sediments
206
 
    if (have_option("/material_phase[0]/sediment")) then
207
 
        call sediment_init()
208
 
    end if
209
 
 
210
205
    ! Initialise Hyperlight
211
206
#ifdef HAVE_HYPERLIGHT
212
207
    if (have_option("ocean_biology/lagrangian_ensemble/hyperlight")) then
343
338
          ! Initialise the OriginalDistanceToBottom field used for wetting and drying
344
339
          if (have_option("/mesh_adaptivity/mesh_movement/free_surface/wetting_and_drying")) then
345
340
             call insert_original_distance_to_bottom(state(1))
346
 
             ! Wetting and drying only works with no poisson guess ... lets check that
347
 
             call get_option("/material_phase::water/scalar_field::Pressure/prognostic/scheme/poisson_pressure_solution", option_buffer)
 
341
             ! Wetting and drying only works with no poisson guess ... let's check that
 
342
             call get_option("/material_phase[0]/scalar_field::Pressure/prognostic/scheme/poisson_pressure_solution", option_buffer)
348
343
             if (.not. trim(option_buffer) == "never") then 
349
 
               FLExit("Please choose 'never' under /material_phase::water/scalar_field::Pressure/prognostic/scheme/poisson_pressure_solution when using wetting and drying")
 
344
               FLExit("Please choose 'never' under /material_phase[0]/scalar_field::Pressure/prognostic/scheme/poisson_pressure_solution when using wetting and drying")
350
345
             end if
351
346
          end if
352
347
       end if
463
458
        call gls_init(state(1))
464
459
    end if
465
460
 
466
 
    ! Initialise k_epsilon
467
 
    have_k_epsilon = .false.
468
 
    keps_option_path="/material_phase[0]/subgridscale_parameterisations/k-epsilon/"
469
 
    if (have_option(trim(keps_option_path))) then
470
 
        have_k_epsilon = .true.
471
 
        call keps_init(state(1))
472
 
    end if
473
 
 
474
461
    ! ******************************
475
462
    ! *** Start of timestep loop ***
476
463
    ! ******************************
642
629
             call solids(state(1), its, nonlinear_iterations)
643
630
          end if
644
631
 
 
632
          ! Do we have the k-epsilon turbulence model?
 
633
          ! If we do then we want to calculate source terms and diffusivity for the k and epsilon 
 
634
          ! fields and also tracer field diffusivities at n + theta_nl
 
635
          do i= 1, size(state)
 
636
             if(have_option("/material_phase["//&
 
637
                  int2str(i-1)//"]/subgridscale_parameterisations/k-epsilon")) then
 
638
                if(timestep == 1 .and. its == 1 .and. have_option('/physical_parameters/gravity')) then
 
639
                   ! The very first time k-epsilon is called, VelocityBuoyancyDensity
 
640
                   ! is set to zero until calculate_densities is called in the momentum equation
 
641
                   ! solve. Calling calculate_densities here is a work-around for this problem.  
 
642
                   sfield => extract_scalar_field(state, 'VelocityBuoyancyDensity')
 
643
                   if(option_count("/material_phase/vector_field::Velocity/prognostic") > 1) then 
 
644
                      call get_phase_submaterials(state, i, submaterials)
 
645
                      call calculate_densities(submaterials, buoyancy_density=sfield)
 
646
                      deallocate(submaterials)
 
647
                   else
 
648
                      call calculate_densities(state, buoyancy_density=sfield)
 
649
                   end if
 
650
                   ewrite_minmax(sfield)
 
651
                end if
 
652
                call keps_advdif_diagnostics(state(i))
 
653
             end if
 
654
          end do
 
655
 
645
656
          field_loop: do it = 1, ntsol
646
657
             ewrite(2, "(a,i0,a,i0)") "Considering scalar field ", it, " of ", ntsol
647
658
             ewrite(1, *) "Considering scalar field " // trim(field_name_list(it)) // " in state " // trim(state(field_state_list(it))%name)
655
666
                end if
656
667
             end if
657
668
 
658
 
             ! do we have the k-epsilon 2 equation turbulence model?
659
 
             if(have_k_epsilon .and. have_option(trim(keps_option_path)//"/scalar_field::"//trim(field_name_list(it)//"/prognostic"))) then
660
 
                if( (trim(field_name_list(it))=="TurbulentKineticEnergy")) then
661
 
                    call keps_tke(state(1))
662
 
                else if( (trim(field_name_list(it))=="TurbulentDissipation")) then
663
 
                    call keps_eps(state(1))
664
 
                endif
665
 
             end if
666
 
 
667
 
            ! Calculate the meltrate
668
 
            if(have_option("/ocean_forcing/iceshelf_meltrate/Holland08/") ) then
 
669
             ! Calculate the meltrate
 
670
             if(have_option("/ocean_forcing/iceshelf_meltrate/Holland08/") ) then
669
671
                if( (trim(field_name_list(it))=="MeltRate")) then
670
 
                    call melt_surf_calc(state(1))
 
672
                   call melt_surf_calc(state(1))
671
673
                endif
672
 
            end if
673
 
 
 
674
             end if
674
675
 
675
676
             call get_option(trim(field_optionpath_list(it))//&
676
677
                  '/prognostic/equation[0]/name', &
677
678
                  option_buffer, default="UnknownEquationType")
678
679
             select case(trim(option_buffer))
679
 
             case ( "AdvectionDiffusion", "ConservationOfMass", "ReducedConservationOfMass", "InternalEnergy", "HeatTransfer" )
 
680
             case ( "AdvectionDiffusion", "ConservationOfMass", "ReducedConservationOfMass", "InternalEnergy", "HeatTransfer", "KEpsilon" )
680
681
                use_advdif=.true.
681
682
             case default
682
683
                use_advdif=.false.
719
720
                else if(have_option(trim(field_optionpath_list(it)) // &
720
721
                     & "/prognostic/spatial_discretisation/continuous_galerkin")) then
721
722
 
722
 
                   call solve_field_equation_cg(field_name_list(it), state(field_state_list(it)), dt)
 
723
                   call solve_field_equation_cg(field_name_list(it), state, field_state_list(it), dt)
723
724
                else
724
725
 
725
726
                   ewrite(2, *) "Not solving scalar field " // trim(field_name_list(it)) // " in state " // trim(state(field_state_list(it))%name) //" in an advdif-like subroutine."
736
737
          if( have_option("/material_phase[0]/subgridscale_parameterisations/GLS/option")) then
737
738
            call gls_diffusivity(state(1))
738
739
          end if
739
 
 
740
 
          ! k_epsilon after the solve on Epsilon has finished
741
 
          if(have_k_epsilon .and. have_option(trim(keps_option_path)//"/scalar_field::ScalarEddyViscosity/diagnostic")) then
742
 
            ! Update the diffusivity, at each iteration.
743
 
            call keps_eddyvisc(state(1))
744
 
          end if
745
740
          
746
741
          !BC for ice melt
747
742
          if (have_option('/ocean_forcing/iceshelf_meltrate/Holland08/calculate_boundaries')) then
929
924
       end if
930
925
 
931
926
    end do timestep_loop
 
927
 
932
928
    ! ****************************
933
929
    ! *** END OF TIMESTEP LOOP ***
934
930
    ! ****************************
947
943
        call gls_cleanup()
948
944
    end if
949
945
 
950
 
    ! cleanup k_epsilon
951
 
    if (have_k_epsilon) then
952
 
        call keps_cleanup()
953
 
    end if
954
 
 
955
 
    if (have_option("/material_phase[0]/sediment")) then
956
 
        call sediment_cleanup()
957
 
    end if
958
 
 
959
946
    ! closing .stat, .convergence and .detector files
960
947
    call close_diagnostic_files()
961
948
 
1033
1020
        call gls_cleanup() ! deallocate everything
1034
1021
    end if
1035
1022
 
1036
 
    ! k_epsilon - we need to deallocate all module-level fields or the memory
1037
 
    ! management system complains
1038
 
    if (have_option("/material_phase[0]/subgridscale_parameterisations/k-epsilon/")) then
1039
 
        call keps_cleanup() ! deallocate everything
1040
 
    end if
1041
 
 
1042
1023
    ! deallocate sub-state
1043
1024
    if(use_sub_state()) then
1044
1025
      do ss = 1, size(sub_state)
1055
1036
    real, intent(inout) :: dt
1056
1037
    integer, intent(inout) :: nonlinear_iterations, nonlinear_iterations_adapt
1057
1038
    type(state_type), dimension(:), pointer :: sub_state
1058
 
    character(len=OPTION_PATH_LEN) :: keps_option_path
 
1039
    integer :: i
1059
1040
 
1060
1041
    ! Overwrite the number of nonlinear iterations if the option is switched on
1061
1042
    if(have_option("/timestepping/nonlinear_iterations/nonlinear_iterations_at_adapt")) then
1120
1101
        call gls_adapt_mesh(state(1))
1121
1102
    end if
1122
1103
 
1123
 
    ! k_epsilon
1124
 
    keps_option_path="/material_phase[0]/subgridscale_parameterisations/k-epsilon/"
1125
 
    if (have_option(trim(keps_option_path)//"/scalar_field::TurbulentKineticEnergy/prognostic") &
1126
 
        &.and. have_option(trim(keps_option_path)//"/scalar_field::TurbulentDissipation/prognostic") &
1127
 
        &.and. have_option(trim(keps_option_path)//"/scalar_field::ScalarEddyViscosity/diagnostic")) then
1128
 
        call keps_adapt_mesh(state(1))
1129
 
    end if
1130
 
 
1131
1104
  end subroutine update_state_post_adapt
1132
1105
 
1133
1106
  subroutine fluids_module_check_options
1219
1192
  end subroutine check_old_code_path
1220
1193
 
1221
1194
  end module fluids_module
1222