diff --git a/src/float_vgrid.f90 b/src/float_vgrid.f90 index f4dd3b9..793d151 100644 --- a/src/float_vgrid.f90 +++ b/src/float_vgrid.f90 @@ -18,12 +18,12 @@ program float_vgrid ' --vgrid vertical grid (default ''ocean_vgrid.nc'') ', & ''] - call set_args('--vgrid "ocean_vgrid.nc"', help_text) + call set_args('--vgrid "ocean_vgrid.nc" --vgrid_type "mom5"', help_text) file = sget('vgrid') call check_file_exist(file) - vgrid = vgrid_t(file) + vgrid = vgrid_t(file, sget('vgrid_type')) call vgrid%float() call vgrid%update_history(get_mycommand()) call vgrid%write(file) diff --git a/src/topography.f90 b/src/topography.f90 index 4545677..c562b57 100644 --- a/src/topography.f90 +++ b/src/topography.f90 @@ -434,9 +434,9 @@ subroutine topography_deseas(this) end subroutine topography_deseas !------------------------------------------------------------------------- - subroutine topography_min_max_depth(this, vgrid_file, level) + subroutine topography_min_max_depth(this, vgrid_file, vgrid_type, level) class(topography_t), intent(inout) :: this - character(len=*), intent(in) :: vgrid_file + character(len=*), intent(in) :: vgrid_file, vgrid_type integer, intent(in) :: level integer(int32) :: i,j @@ -450,7 +450,7 @@ subroutine topography_min_max_depth(this, vgrid_file, level) this%min_level = level - vgrid = vgrid_t(vgrid_file) + vgrid = vgrid_t(vgrid_file, vgrid_type) this%min_depth = vgrid%zeta(this%min_level) this%max_depth = vgrid%zeta(vgrid%nlevels) @@ -498,9 +498,9 @@ subroutine topography_fill_fraction(this, sea_area_fraction) end subroutine topography_fill_fraction !------------------------------------------------------------------------- - subroutine topography_nonadvective(this, vgrid_file, potholes, coastal, fix) + subroutine topography_nonadvective(this, vgrid_file, vgrid_type, potholes, coastal, fix) class(topography_t), intent(inout) :: this - character(len=*), intent(in) :: vgrid_file + character(len=*), intent(in) :: vgrid_file, vgrid_type logical, intent(in) :: potholes, coastal, fix real(real32), allocatable :: depth_halo(:,:) @@ -516,7 +516,7 @@ subroutine topography_nonadvective(this, vgrid_file, potholes, coastal, fix) integer(int32) :: im, ip, jm, jp integer(int32) :: nseas - vgrid = vgrid_t(vgrid_file) + vgrid = vgrid_t(vgrid_file, vgrid_type) write(output_unit,*) 'Zeta dimensions', 2*vgrid%nlevels + 1, vgrid%nlevels allocate(zw(0:vgrid%nlevels)) zw = real(vgrid%zeta) diff --git a/src/topogtools.f90 b/src/topogtools.f90 index de8c075..3bde804 100644 --- a/src/topogtools.f90 +++ b/src/topogtools.f90 @@ -11,7 +11,7 @@ program topogtools character(len=:), allocatable :: help_general(:), help_gen_topo(:), help_deseas(:), help_min_max_depth(:) character(len=:), allocatable :: help_fill_fraction(:), help_fix_nonadvective(:), help_check_nonadvective(:), help_mask(:) character(len=80) :: version_text(1) - character(len=:), allocatable :: file_in, file_out, hgrid, vgrid, grid_type + character(len=:), allocatable :: file_in, file_out, hgrid, vgrid type(topography_t) :: topog real(real32) :: sea_area_fraction integer :: ii @@ -55,14 +55,16 @@ program topogtools ''] help_min_max_depth = [character(len=80) :: & 'usage: topogtools min_max_depth --input --output ', & - ' --level [--vgrid ] ', & + ' --level ', & + ' [--vgrid --vgrid_type ] ', & ' ', & 'Set minimum depth to the depth at a specified level and set maximum depth to ', & 'deepest in . is the minimum number of depth levels (e.g. 4). ', & 'Can produce non-advective cells. ', & ' ', & 'Options ', & - ' --vgrid vertical grid (default ''ocean_vgrid.nc'') ', & + ' --vgrid vertical grid (default ''ocean_vgrid.nc'') ', & + ' --vgrid_type can be ''mom5'' or ''mom6'' (default ''mom5'') ', & ''] help_fill_fraction = [character(len=80) :: & 'usage: topogtools fill_fraction --input --output ', & @@ -81,6 +83,7 @@ program topogtools ' ', & 'Options ', & ' --vgrid vertical grid (default ''ocean_vgrid.nc'') ', & + ' --vgrid_type can be ''mom5'' or ''mom6'' (default ''mom5'') ', & ' --potholes fix potholes ', & ' --coastal-cells fix non-advective coastal cells ', & ''] @@ -94,6 +97,7 @@ program topogtools ' ', & 'Options ', & ' --vgrid vertical grid (default ''ocean_vgrid.nc'') ', & + ' --vgrid_type can be ''mom5'' or ''mom6'' (default ''mom5'') ', & ' --potholes check for potholes ', & ' --coastal-cells check for non-advective coastal cells ', & ''] @@ -112,14 +116,15 @@ program topogtools case ('deseas') call set_args('--input:i "unset" --output:o "unset"', help_deseas, version_text) case ('min_max_depth') - call set_args('--input:i "unset" --output:o "unset" --vgrid "ocean_vgrid.nc" --level 0', help_min_max_depth, version_text) + call set_args('--input:i "unset" --output:o "unset" --vgrid "ocean_vgrid.nc" --vgrid_type "mom5" --level 0', & + help_min_max_depth, version_text) case('fill_fraction') call set_args('--input:i "unset" --output:o "unset" --fraction 0.0', help_fill_fraction, version_text) case ('fix_nonadvective') - call set_args('--input:i "unset" --output:o "unset" --vgrid "ocean_vgrid.nc" --potholes F --coastal-cells F', & - help_fix_nonadvective, version_text) + call set_args('--input:i "unset" --output:o "unset" --vgrid "ocean_vgrid.nc" --vgrid_type "mom5" --potholes F & + &--coastal-cells F', help_fix_nonadvective, version_text) case ('check_nonadvective') - call set_args('--input:i "unset" --vgrid "ocean_vgrid.nc" --potholes F --coastal-cells F', & + call set_args('--input:i "unset" --vgrid "ocean_vgrid.nc" --vgrid_type "mom5" --potholes F --coastal-cells F', & help_check_nonadvective, version_text) case ('mask') call set_args('--input:i "unset" --output:o "unset"', help_mask, version_text) @@ -174,7 +179,7 @@ program topogtools case ('min_max_depth') topog = topography_t(file_in) - call topog%min_max_depth(vgrid, iget('level')) + call topog%min_max_depth(vgrid, sget('vgrid_type'), iget('level')) call topog%update_history(get_mycommand()) call topog%write(file_out) @@ -191,13 +196,13 @@ program topogtools case ('fix_nonadvective') topog = topography_t(file_in) - call topog%nonadvective(vgrid, lget('potholes'), lget('coastal-cells'), fix=.true.) + call topog%nonadvective(vgrid, sget('vgrid_type'), lget('potholes'), lget('coastal-cells'), fix=.true.) call topog%update_history(get_mycommand()) call topog%write(file_out) case ('check_nonadvective') topog = topography_t(file_in) - call topog%nonadvective(vgrid, lget('potholes'), lget('coastal-cells'), fix=.false.) + call topog%nonadvective(vgrid, sget('vgrid_type'), lget('potholes'), lget('coastal-cells'), fix=.false.) case ('mask') topog = topography_t(file_in) diff --git a/src/vgrid.f90 b/src/vgrid.f90 index a33d75a..875d25c 100644 --- a/src/vgrid.f90 +++ b/src/vgrid.f90 @@ -9,6 +9,7 @@ module vgrid type vgrid_t ! Vertical grid variable + character(len=:), allocatable :: type integer :: nlevels = 0 real(real64), allocatable :: zeta(:) real(real64), allocatable :: zeta_super(:) @@ -31,14 +32,22 @@ module vgrid contains !------------------------------------------------------------------------- - type(vgrid_t) function vgrid_constructor(filename) result(vgrid) - character(len=*), intent(in) :: filename + type(vgrid_t) function vgrid_constructor(filename, type) result(vgrid) + character(len=*), intent(in) :: filename, type integer(int32) :: ncid, zeta_id, did(1), zeta_len, author_len, history_len ! NetCDF ids and dims + real(real64), allocatable :: zeta_var(:) write(output_unit,'(3a)') "Reading vgrid from file '", trim(filename), "'" vgrid%original_file = filename + vgrid%type = type + select case (type) + case ("mom5", "mom6") + case default + write(error_unit,'(3a)') "ERROR: '", type, "' is not a valid vertical grid type." + error stop + end select ! Open file call handle_error(nf90_open(trim(filename), nf90_nowrite, ncid)) @@ -46,14 +55,11 @@ type(vgrid_t) function vgrid_constructor(filename) result(vgrid) ! Get dimension call handle_error(nf90_inq_dimid(ncid, 'nzv', did(1))) call handle_error(nf90_inquire_dimension(ncid, did(1), len=zeta_len)) - vgrid%nlevels = zeta_len/2 ! Get zeta - allocate(vgrid%zeta_super(zeta_len)) - allocate(vgrid%zeta(0:vgrid%nlevels)) + allocate(zeta_var(zeta_len)) call handle_error(nf90_inq_varid(ncid, 'zeta', zeta_id)) - call handle_error(nf90_get_var(ncid, zeta_id, vgrid%zeta_super)) - vgrid%zeta = vgrid%zeta_super(1:zeta_len:2) + call handle_error(nf90_get_var(ncid, zeta_id, zeta_var)) if (nf90_inquire_attribute(ncid, zeta_id, 'author', len=author_len) == nf90_noerr) then allocate(character(len=author_len) :: vgrid%author) call handle_error(nf90_get_att(ncid, zeta_id, 'author', vgrid%author)) @@ -65,6 +71,21 @@ type(vgrid_t) function vgrid_constructor(filename) result(vgrid) call handle_error(nf90_get_att(ncid, nf90_global, 'history', vgrid%history)) end if + ! Handle the different types of grids + select case (type) + case ("mom5") + vgrid%nlevels = zeta_len/2 + allocate(vgrid%zeta_super(zeta_len)) + allocate(vgrid%zeta(0:vgrid%nlevels)) + vgrid%zeta_super = zeta_var + vgrid%zeta = zeta_var(1:zeta_len:2) + + case ("mom6") + vgrid%nlevels = zeta_len - 1 + allocate(vgrid%zeta(0:vgrid%nlevels)) + vgrid%zeta = zeta_var + end select + ! Close file call handle_error(nf90_close(ncid)) @@ -151,13 +172,15 @@ subroutine vgrid_float(this) real(real32), allocatable :: zeta_float(:) - allocate(zeta_float(2*this%nlevels+1)) - zeta_float = this%zeta_super - this%zeta_super = zeta_float - deallocate(zeta_float) + if (allocated(this%zeta_super)) then + allocate(zeta_float(2*this%nlevels+1)) + zeta_float = real(this%zeta_super, real32) + this%zeta_super = zeta_float + deallocate(zeta_float) + end if allocate(zeta_float(0:this%nlevels)) - zeta_float = this%zeta + zeta_float = real(this%zeta, real32) this%zeta = zeta_float deallocate(zeta_float)