Skip to content

Commit

Permalink
Disabled reading of grid information, since current version is outdat…
Browse files Browse the repository at this point in the history 
…ed and not consistent with input files. Removed unused parameter ydyn.beta_gl_sep. Reduced initial timestep for SIA solver for 0.1yr equilibrium, before moving forward. Also updated ncio/nml to latest versions (though not really necessary in the end).
  • Loading branch information
Alexander Robinson committed Dec 18, 2021
1 parent bd5ed43 commit 33bd3e1
Show file tree
Hide file tree
Showing 7 changed files with 185 additions and 36 deletions.
180 changes: 155 additions & 25 deletions libs/ncio.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@

module ncio

use ieee_arithmetic
use netcdf

implicit none
Expand All @@ -36,6 +37,8 @@ module ncio

character(len=NC_STRLEN), parameter :: NC_CHARDIM = "strlen"

character(len=3), parameter :: GRID_MAPPING_NAME_DEFAULT = "crs"

type ncvar
character (len=NC_STRLEN) :: name, long_name, standard_name, units
character (len=NC_STRLEN) :: axis, calendar, grid_mapping
Expand Down Expand Up @@ -381,6 +384,8 @@ subroutine nc4_read_internal_numeric(filename,name,dat,size_in,start,count,xtype
integer, optional :: iostat
integer :: i, status

double precision, parameter :: missing_value_default = -9999.0

! Open the file.
call nc_check_open(filename, ncid, nf90_nowrite, nc_id)

Expand Down Expand Up @@ -440,6 +445,15 @@ subroutine nc4_read_internal_numeric(filename,name,dat,size_in,start,count,xtype
if (.not. present(ncid)) call nc_check( nf90_close(nc_id) )

if (v%missing_set) then

! === SPECIAL CASE: missing_value == NaN ====

! Replace NaNs with internal missing value to avoid crashes
where ( ieee_is_nan(dat) ) dat = missing_value_default
v%missing_value = missing_value_default

! ===========================================

where( dabs(dat-v%missing_value) .gt. NC_TOL ) dat = dat*v%scale_factor + v%add_offset

! Fill with user-desired missing value
Expand Down Expand Up @@ -543,12 +557,12 @@ subroutine nc_v_init(v,name,xtype,ndims_in,long_name,standard_name, &
v%missing_set = .TRUE.
v%missing_value = -9999d0
v%FillValue = v%missing_value
v%FillValue_set = .FALSE.
v%FillValue_set = .TRUE.

v%xtype = "NF90_DOUBLE"
v%coord = .FALSE.

v%grid_mapping = ""
v%grid_mapping = GRID_MAPPING_NAME_DEFAULT

! If args are present, reassign these variables
if ( present(long_name) ) v%long_name = trim(long_name)
Expand Down Expand Up @@ -880,14 +894,15 @@ subroutine nc_put_att(ncid, v)
if (v%FillValue_set) then
select case(trim(v%xtype))
case("NF90_INT")
call nc_check( nf90_put_att(ncid, v%varid, "FillValue", int(v%FillValue)) )
call nc_check( nf90_put_att(ncid, v%varid, "_FillValue", int(v%FillValue)) )
case("NF90_FLOAT")
call nc_check( nf90_put_att(ncid, v%varid, "FillValue", real(v%FillValue)) )
call nc_check( nf90_put_att(ncid, v%varid, "_FillValue", real(v%FillValue)) )
case("NF90_DOUBLE")
call nc_check( nf90_put_att(ncid, v%varid, "FillValue", v%FillValue) )
call nc_check( nf90_put_att(ncid, v%varid, "_FillValue", v%FillValue) )
end select
end if


end if

! Add additional variable attributes
Expand Down Expand Up @@ -1014,10 +1029,17 @@ subroutine nc_get_att(ncid, v, readmeta)
call nc_get_att_double(ncid,v%varid,"missing_value",v%missing_value,stat)
if (stat .eq. noerr) v%missing_set = .TRUE.

stat = nc_check_att( nf90_get_att(ncid, v%varid, "FillValue", tmpi) )
stat = nc_check_att( nf90_get_att(ncid, v%varid, "_FillValue", tmpi) )
if (stat .eq. noerr) then
v%FillValue = dble(tmpi)
v%FillValue_set = .TRUE.

! ajr, 2020-08-20
! Overwrite missing value with Fillvalue, since missing_value
! attribute has been deprecated.
v%missing_value = v%FillValue
v%missing_set = .TRUE.

end if

case DEFAULT
Expand All @@ -1034,7 +1056,7 @@ subroutine nc_get_att(ncid, v, readmeta)
call nc_get_att_double(ncid,v%varid,"missing_value",v%missing_value,stat)
if (stat .eq. noerr) v%missing_set = .TRUE.

stat = nc_check_att( nf90_get_att(ncid, v%varid, "FillValue", tmp) )
stat = nc_check_att( nf90_get_att(ncid, v%varid, "_FillValue", tmp) )
if (stat .eq. noerr) then
v%FillValue = tmp
v%FillValue_set = .TRUE.
Expand Down Expand Up @@ -1418,61 +1440,166 @@ subroutine nc_write_attr_global_dp(filename, name, value)
call nc_check( nf90_close(ncid) )
end subroutine

subroutine nc_write_map(filename,name,lambda,phi,alpha,x_e,y_n, ncid)
subroutine nc_write_map(filename,grid_mapping_name,lambda,phi,alpha,x_e,y_n, &
is_sphere,semi_major_axis,inverse_flattening,ncid)
! CF map conventions can be found here:
! http://cfconventions.org/Data/cf-conventions/cf-conventions-1.6/build/cf-conventions.html#appendix-grid-mappings

implicit none

character(len=*) :: filename, name

character(len=*), intent(IN) :: filename, grid_mapping_name
double precision, intent(IN), optional :: lambda, phi
double precision, intent(IN), optional :: alpha, x_e, y_n
logical, intent(IN), optional :: is_sphere
double precision, intent(IN), optional :: semi_major_axis
double precision, intent(IN), optional :: inverse_flattening
integer, intent(in), optional :: ncid

! Local variables
integer :: nc_id, varid, stat
integer, intent(in), optional :: ncid
double precision, optional :: lambda, phi, alpha, x_e, y_n
character(len=56) :: crs_name
double precision :: phi_proj_orig

integer, parameter :: noerr = NF90_NOERR

! Define general name of integer variable to
! hold coordinate reference system (CRS) information
! Always, use crs by default to be consistent with others.
crs_name = "crs"

! Open the file, set for redefinition
call nc_check_open(filename, ncid, nf90_write, nc_id)
call nc_check( nf90_redef(nc_id) )

! Check if grid mapping has been defined in this file
! (if not, define it according to input arguments)
stat = nf90_inq_varid(nc_id, trim(name), varid)
stat = nf90_inq_varid(nc_id, trim(crs_name), varid)

if ( stat .ne. noerr ) then
! Define the mapping variable as an integer with no dimensions,
! and include the grid mapping name
call nc_check( nf90_def_var(nc_id, trim(name), NF90_INT, varid) )
call nc_check( nf90_put_att(nc_id,varid, "grid_mapping_name", trim(name)) )


call nc_check( nf90_def_var(nc_id, trim(crs_name), NF90_INT, varid) )
! Add grid attributes depending on grid_mapping type
select case(trim(name))
select case(trim(grid_mapping_name))

case("stereographic")
! Including 'oblique_stereographic'

if ( (.not. present(lambda)) .or. &
(.not. present(phi)) .or. &
(.not. present(alpha)) .or. &
(.not. present(x_e)) .or. &
(.not. present(y_n)) ) then

write(*,"(a,a)") "ncio:: nc_write_map:: Error: ", &
"All grid_mapping arguments must be provided ", &
"(lambda, phi, alpha, x_e, y_n)."
stop

end if

! Add grid mapping attributes
call nc_check( nf90_put_att(nc_id,varid, "grid_mapping_name", trim(grid_mapping_name)) )
call nc_check( nf90_put_att(nc_id,varid, "longitude_of_projection_origin", lambda) )
call nc_check( nf90_put_att(nc_id,varid, "latitude_of_projection_origin", phi) )
if (present(alpha)) &
call nc_check( nf90_put_att(nc_id,varid, "angle_of_oblique_tangent", alpha) )
call nc_check( nf90_put_att(nc_id,varid, "scale_factor_at_projection_origin", 1.d0) )
call nc_check( nf90_put_att(nc_id,varid, "false_easting", x_e) )
call nc_check( nf90_put_att(nc_id,varid, "false_northing", y_n) )

case("polar_stereographic")

if ( (.not. present(lambda)) .or. &
(.not. present(phi)) .or. &
(.not. present(x_e)) .or. &
(.not. present(y_n)) ) then

write(*,"(a,a)") "ncio:: nc_write_map:: Error: ", &
"All grid_mapping arguments must be provided ", &
"(lambda, phi, x_e, y_n)."
stop

end if

! Determine latitude_of_projection_origin, since it must
! be either -90 or +90 for a polar_stereographic projection:
if (phi .gt. 0.0d0) then
phi_proj_orig = 90.0d0
else
phi_proj_orig = -90.0d0
end if

! Add grid mapping attributes
call nc_check( nf90_put_att(nc_id,varid, "grid_mapping_name", trim(grid_mapping_name)) )
call nc_check( nf90_put_att(nc_id,varid, "straight_vertical_longitude_from_pole", lambda) )
call nc_check( nf90_put_att(nc_id,varid, "latitude_of_projection_origin", phi) )
if (present(alpha)) &
call nc_check( nf90_put_att(nc_id,varid, "angle_of_oblique_tangent", alpha) )
call nc_check( nf90_put_att(nc_id,varid, "scale_factor_at_projection_origin", 1.d0) )
call nc_check( nf90_put_att(nc_id,varid, "latitude_of_projection_origin", phi_proj_orig) )
call nc_check( nf90_put_att(nc_id,varid, "standard_parallel", phi) )

call nc_check( nf90_put_att(nc_id,varid, "false_easting", x_e) )
call nc_check( nf90_put_att(nc_id,varid, "false_northing", y_n) )

case("lambert_conformal_conic")

call nc_check( nf90_put_att(nc_id,varid, "longitude_of_central_meridian", lambda) )
call nc_check( nf90_put_att(nc_id,varid, "latitude_of_projection_origin", phi) )
if (present(alpha)) &
call nc_check( nf90_put_att(nc_id,varid, "standard_parallel", alpha) )

case("latitude_longitude","latlon","gaussian")

! Pass - no grid mapping needed for a latitude_longitude grid

! Add grid mapping attributes
!call nc_check( nf90_put_att(nc_id,varid, "grid_mapping_name", "latitude_longitude") )
! == No additional parameters needed ==

case DEFAULT
! Do nothing

end select
end select


select case(trim(grid_mapping_name))

case("latitude_longitude","latlon","gaussian")

write(*,"(a,a,a)") "ncio:: nc_write_map:: ",trim(filename)//" : ",trim(name)
! Pass - do nothing for latitude_longitude grids

case DEFAULT
! For other projections, add planet information if available

! Add planet information if desired
if (present(is_sphere) .and. &
present(semi_major_axis) .and. &
present(inverse_flattening)) then

if (is_sphere) then
call nc_check( nf90_put_att(nc_id,varid, "semi_major_axis", semi_major_axis) )
call nc_check( nf90_put_att(nc_id,varid, "inverse_flattening", 0.0d0) )
else
call nc_check( nf90_put_att(nc_id,varid, "semi_major_axis", semi_major_axis) )
call nc_check( nf90_put_att(nc_id,varid, "inverse_flattening", inverse_flattening) )

end if

else if (present(is_sphere) .or. &
present(semi_major_axis) .or. &
present(inverse_flattening)) then

write(*,*) "ncio:: nc_write_map:: Error: to write planet information, &
&all three planet parameters must be provided: &
&is_sphere, semi_major_axis, inverse_flattening. &
&Try again."

stop
end if

end select

write(*,"(a,a,a)") "ncio:: nc_write_map:: ",trim(filename)//" : ",trim(grid_mapping_name)

end if

Expand All @@ -1495,7 +1622,6 @@ end subroutine nc_write_map
!! @param units NetCDF attribute of the units of the variable (optional)
!! @param axis NetCDF attribute of the standard axis of the variable (optional)
!! @param calendar NetCDF attribute of the calendar type to be used for time dimensions (optional)

subroutine nc_write_dim_int_pt(filename,name,x,dx,nx, &
long_name,standard_name,units,axis,calendar,unlimited, ncid)

Expand Down Expand Up @@ -1723,6 +1849,10 @@ subroutine nc_write_dim_internal(filename,name,xtype,x, &

call nc_v_init(v,name=trim(name),xtype=xtype,coord=.TRUE.)

! Ensure grid_mapping is set to blank since a dimension cannot
! have a grid mapping.
v%grid_mapping = ""

!! Now fill in values of arguments that are present
if ( present(long_name) ) v%long_name = trim(long_name)
if ( present(standard_name) ) v%standard_name = trim(standard_name)
Expand Down
29 changes: 23 additions & 6 deletions libs/nml.f90
View file
Original file line number Diff line number Diff line change
@@ -1,10 +1,14 @@

module nml

use, intrinsic :: iso_fortran_env, only : input_unit, output_unit, error_unit

implicit none

logical :: VERBOSE = .FALSE. !! should freshly read namelist be printed to screen?
logical :: VERBOSE = .FALSE. !! should freshly read namelist be printed to screen?
logical :: ERROR_NO_PARAM = .TRUE. !! Should error be thrown if parameter isn't found?

integer, parameter :: io_unit_err = error_unit

interface nml_read
module procedure nml_read_string, nml_read_double, nml_read_float
Expand Down Expand Up @@ -46,7 +50,7 @@ subroutine nml_set_verbose(verb)

return

end subroutine nml_set_verbose
end subroutine nml_set_verbose

subroutine nml_replace(s,text,rep,outs)
! Adapted from FUNCTION Replace_Text:
Expand Down Expand Up @@ -90,10 +94,12 @@ subroutine nml_read_internal(filename,group,name,value,comment)
character(len=*), intent(IN) :: filename, group, name
character(len=*), intent(INOUT), optional :: comment

! Local variables
integer :: io, file_unit
integer :: iostat, l, ltype
character(len=1000) :: line, name1, value1, comment1
logical :: ingroup
logical :: found_param

! Open the nml filename to be read, or get file units if already open
inquire(file=trim(filename),NUMBER=file_unit)
Expand All @@ -111,7 +117,8 @@ subroutine nml_read_internal(filename,group,name,value,comment)
end if


ingroup = .FALSE.
ingroup = .FALSE.
found_param = .FALSE.

do l = 1, 5000
read(io,"(a1000)",iostat=iostat) line
Expand All @@ -121,8 +128,9 @@ subroutine nml_read_internal(filename,group,name,value,comment)
! Check if the parameter has been found
if (ingroup .and. ltype == 3) then
if (trim(name1) == trim(name)) then
value = trim(value1)
comment = trim(comment1)
value = trim(value1)
comment = trim(comment1)
found_param = .TRUE.
exit
end if
end if
Expand All @@ -141,6 +149,15 @@ subroutine nml_read_internal(filename,group,name,value,comment)
close(io)
endif

if (ERROR_NO_PARAM .and. (.not. found_param)) then
write(io_unit_err,*) ""
write(io_unit_err,*) "nml:: Error: parameter not found."
write(io_unit_err,*) "Filename: ", trim(filename)
write(io_unit_err,*) "Group: ", trim(group)
write(io_unit_err,*) "Parameter: ", trim(name)
stop "Program stopped."
end if

return

end subroutine nml_read_internal
Expand Down Expand Up @@ -170,7 +187,7 @@ subroutine nml_read_string(filename,group,name,value,comment,init)

return

end subroutine nml_read_string
end subroutine nml_read_string

subroutine nml_read_double(filename,group,name,value,comment,init)

Expand Down
1 change: 1 addition & 0 deletions par/yelmo_initmip.nml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -94,6 +94,7 @@
beta_gl_stag = 1 ! 0: simple staggering, 1: Upstream beta at gl, 2: downstream, 3: f_grnd_ac scaling
beta_gl_f = 1.0 ! [-] Scaling of beta at the grounding line (for beta_gl_scale=0)
taud_gl_method = 0 ! 0: binary, no subgrid, 1: Two-sided gradient
calc_diffusivity = False
H_grnd_lim = 500.0 ! [m] For beta_gl_scale=1, reduce beta linearly between H_grnd=0 and H_grnd_lim
H_sed_sat = 250.0 ! [m] Sediment thickness at which sediment effect is saturated
cb_method = 1 ! -1: set externally; 1: calculate cb online
Expand Down
Loading

0 comments on commit 33bd3e1

Please sign in to comment.